Univerza v Ljubljani Fakulteta za elektrotehniko
Blaž Valič
Problem velikih geometrijskih
razmerij pri modeliranju bioloških
tkiv z metodo končnih elementov
Magistrska naloga
Mentor: prof. dr. Damijan Miklavčič Ljubljana, junij 2004
Zahvala
Gotovo ima veliko zaslug za nastanek te naloge prof. dr. Damijan Miklavčič, moj mentor tako pri podiplomskem študiju kot tudi pri pripravi magistrske naloge. Vedno me je znal spodbujati in mi svetovati, predvsem takrat, ko se je zataknilo, ko nisem imel več energije oziroma ideje, kako naprej.
Pri gradnji modelov, ki si bili osnova naloge, sem veliko časa presedel za računalnikom in se jezil nad muhami FEMLABA. Pri tem mi je v nekaj mesecih, ki jih je preživel v našem laboratoriju, delal družbo mag. Igor Lacković. Skupaj je bilo lažje, še posebej zato, ker sva hkrati odkrivala lepote in tegobe modelov in algoritmov in                      takšnih       sporočil:
"Out of memory", "Step size too smal", …
A nisva bila sama, takšna sporočila pozna tudi dr. Davorka Šel. Verjamem, da se je tudi ona pri gradnji modelov srečevala s podobnimi težavami, ki mi jih je delno prihranila, saj sem lahko uporabil njeno geometrijo modela glave s tumorjem in elektrodami. Hvala tudi za vse odgovore na moje maile.
V prijetnem okolju je prijetno delati, k temu pa so veliko prispevali sodelavci v Laboratoriju za biokibernetiko. Vsem najlepša hvala za prijetne pogovore, nasvete, ideje, čas in dobro voljo.
Hvala doc. dr. Alenki Maček Lebar za nasvete pri pisanju ter za pregled naloge.
Nalogo je lektorirala prijateljica Sanja Leban. Hvala!
»Last, but not least«, hvala Andreji, Nejcu, Aljažu in Niki, predvsem za potrpljenje in podporo pri mojem delu.
i
Kazalo
KAZALO	I
POVZETEK	III
ABSTRACT	V
1    UVOD	1
1.1    Zakaj modeliranje in numeriČne metode v biomedicini?
1.2    Osnove modeliranja in izraČuna z numeriČnimi metodami
1.2.1    Posebnosti modeliranja v biomedicini
1.3    Problem velikih geometrijskih razmerij
1.4    Možnosti za reševanje problema velikih geometrijskih razmerij
1.5    Cilji naloge                                                                                                                    6
2    TEORIJA                                                                                                                             7
2.1    Elektromagnetno polje
2.1.1    Tokovno polje                                                                                                             10
2.2    Modeliranje in izraČun z numeriČnimi metodami                                                    11
2.2.1    Metoda konČnih diferenc                                                                                           12
2.2.2    Metoda konČnih elementov                                                                                       14
3    REŠEVANJE PROBLEMA VELIKIH GEOMETRIJSKIH RAZMERIJ                           18
3.1    Uporaba skaliranih elektrod                                                                                    18
3.2    Uporaba nelokalnega sklapljanja v FEMLAB-u                                                     20
3.2.1    Model možganov s tumorjem                                                                                     22
3.2.2    Model možganov s tumorjem in nelokalnim sklapljanjem                                      24
3.2.3    Model celice z nelokalnim sklapljanjem                                                                  26
ii
4    REZULTATI                                                                                                                      30
4.1    Uporaba skaliranih elektrod                                                                                    30
4.2    Model možganov s tumorjem                                                                                     35
4.3    Model celice                                                                                                                39
5    RAZPRAVA                                                                                                                       42
5.1    Uporaba skaliranih elektrod                                                                                    42
5.2    Model možganov s tumorjem in model celice                                                         44
6    ZAKLJUČEK                                                                                                                    47
7
LITERATURA
48
iii
Povzetek
Dandanes se modeliranje in izračun z numeričnimi metodami pogosto uporablja v biomedicini in biomedicinski tehniki. Omogoča nam, da lahko na primer preizkušamo občutljivosti sistemov na posamezne parametre, kar je ponavadi veliko dražje, težje ali celo nemogoče eksperimentalno ali analitično, prav tako se lahko pri načrtovanju sistema ali naprave opiramo na podatke, pridobljene iz odzivov modela, in podobno.
V  Laboratoriju za biokibernetiko z modeliranjem in izračun z numeričnimi metodami določamo predvsem gostoto električnega toka in porazdelitev električne poljske jakosti v tkivih in računamo transmembransko napetost na celicah. V obeh primerih se v geometriji modela srečujemo z zelo velikimi geometrijskimi razmerji (razmerje med najmanjšo in največjo dimenzijo v geometriji). To predstavlja programskim paketom za modeliranje in izračun z numeričnimi metodami veliko težavo, saj algoritmi, ki so namenjeni razdelitvi prostora računanja v manjša območja, elemente (delitev prostora imenujemo mrežo), ne uspejo zgraditi mreže v takšni geometriji. Problem smo do sedaj reševali z uporabo nadomestnih prevodnosti za celično suspenzijo, z uporabo simetrij in antisimetrij v modelu ter z omejitvijo področja računanja. Vsak način je imel svoje slabosti ali je bil omejen s področjem uporabe.
V  nalogi smo preučili dve novi možnosti za reševanje problema velikih geometrijskih razmerij, skaliranje elektrod in nelokalno sklapljanje, ki omogoča sklopitev dveh točk, ki geometrijsko nista povezani. V problemih, ki rešujemo v Laboratoriju za biokibernetiko, velikokrat kot vir elektromagnetnega polja nastopajo elektrode, ki imajo relativno majhen premer. Z modeli smo ocenili, kolikokrat lahko povečamo premer elektrod ter temu ustrezno znižamo električni potencial na njih, da bo rezultat še dovolj dober. Z nelokalnim sklapljanjem smo zgradili model možganov s tumorjem in vzporednimi osemigelnimi elektrodami ter model celice z membrano. V obeh modelih smo izrezali tisti del celotne geometrije, ki predstavlja največ težav za gradnjo mreže oziroma smo v njem želeli rešitev z večjo prostorsko ločljivostjo in ga vstavili v drugo, manjšo geometrijo. Robne pogoje na mestu rezanja smo določili tako, da so se preslikali iz večje geometrije na manjšo in obratno.
iv
Ugotovili smo, da je skaliranje elektrod uporabna metoda, ki v primeru dveh igelnih elektrod omogoča povečanje premera elektrod za 2 do 5 krat brez večjih napak v porazdelitvi absolutne vrednosti električne poljske jakosti. V primeru drugačne geometrije elektrod je potrebno tako obliko nadomestnih elektrod kakor tudi znižanje električnega potenciala na nadomestnih elektrodah določiti s pomočjo modelov.
Z nelokalnim sklapljanjem smo uspešno zgradili model možganov s tumorjem in vzporednimi osem igelnimi elektrodami. Sicer je bila rešitev modela mogoča tudi brez uporabe nelokalnega sklapljanja, a je bila gradnja mreže v primeru nelokalnega sklapljanja enostavnejša, kljub temu da nismo iz modela izrezali elektrod, kakor je bilo potrebno v modelu brez nelokalnega sklapljanja, rezultat pa je bil bolj podroben. Z modelom celice z membrano je bilo kljub uporabi nelokalnega sklapljanja več težav, predvsem je slabo konvergiral numerični izračun. Izkazalo se je, da je nelokalno sklapljanje zelo uporabno v nekaterih geometrijah ter omogoča izračune takšnih geometrij, ki bi se jih brez uporabe te metode z nam dosegljivimi programskimi paketi ne dalo.
v
Abstract
Nowadays numerical modeling is a common method in the field of biomedicine and biomedical engineering. It enables to analyze sensitivity of systems on different parameters, which can be harder or even impossible experimentally or analytically; it enables to use data acquired from the model in the design of a system or device.
In the Laboratory of Biocybernetics at the University of Ljubljana, Faculty of Electrical Engineering, we use numerical modeling mainly for calculation of electric current and electric field intensity distribution in tissue, and for calculation of the transmembrane voltage on cells. In both cases there are very large geometric ratios in the geometry of the model (ratio between the smallest and the largest dimension in geometry). Numerical modeling programs often have problems with large geometric ratios, because algorithms for partitioning of the calculating space into smaller areas, elements, are unable to perform partitioning in such geometries. We already used some techniques to avoid such problems. For example, we used effective conductivity for cell suspension, symmetries and anti-symmetries in the models to reduce the size of the models and we also reduced the size of the models by defining isolating boundary conditions in the models. Each of the above technique however has some disadvantages or is limited by the application area.
In this work we analyzed two new techniques, which may help overcoming the problem of large geometric ratio: electrode scaling and non-local coupling, which enables us to connect points, which are not geometrically connected. In the Laboratory of Biocibernetics we often model problems, where electrodes with relatively small diameter are used as electromagnetic field sources. Using electrode scaling we estimated, how many times we can increase the diameter of the electrodes and correspondingly decrease the potential on them, keeping the results good enough. With non-local coupling we build model of brain with parallel eight needles electrodes and model of a cell with the membrane. In both models we cut out a part of the whole geometry, which was most problematic for partitioning or where we wanted better spatial resolution, and inserted it in the second, small geometry in the same model. We
vi
defined coupled boundary conditions on the cut planes between big and small geometry and vice versa.
Our results show that electrode scaling is a useful technique enabling us to increase the diameter of the electrodes by a factor between two and five without major error in the distribution of absolute value of electric field intensity in the case of two needle electrodes. In the case of different electrode geometry it is necessary to determine the shape of scaled electrodes as well as the factor of corresponding decrease of electric potential on the electrodes by new models.
With non-local coupling we successfully built model of brain with parallel eight needles electrodes. It was also possible to build the model without non-local coupling. However partitioning of the space of calculation was much easier if we used non-local coupling as well as we increased the spatial resolution. Despite non-local coupling we had more problems with model of the cell with membrane; mostly the convergence of the numerical calculation was slow. Nevertheless, non-local coupling is very useful in numerous geometries. It enables numerical modeling of geometries, which otherwise can not be solved with available program packages for numerical modeling.
1 Uvod
1
1 Uvod
1.1 Zakaj modeliranje in numerične metode v biomedicini?
Želja po poznavanju procesov v naravi in vplivanju nanje je bila in bo pomembna gonilna sila raziskav in razvoja. Človek se je že v preteklosti in se še danes srečuje s številnimi različnimi procesi v naravi, ki so tako ali drugače pomembni zanj. Tri področja, ki se ukvarjajo s takšnimi procesi in pridobivajo na pomenu v sodobni človeški družbi, so biomedicina, biotehnologija in biomedicinska tehnika. »Področje biomedicina združuje znanja biokemije, molekularne in celične biologije, farmacije, medicine, veterinarske medicine in drugih znanj« [Univerza v Ljubljani, 2003] z namenom izboljševanja kvalitete življenja. »Biotehnologijo v splošnem pomenu opredeljujemo kot integrirano uporabo molekularno-bioloških in inženirskih znanj za komercialne aplikacije organizmov« [Biotehniška fakulteta, Univerza v Ljubljani, 2004]. Biomedicinska tehnika je tehnično področje, ki vključuje znanja biomedicine in klinične prakse. Vključuje: »1. odkrivanje novih znanj in razumevanja živih organizmov s pomočjo inovativnih in samostojnih aplikacij eksperimentalnih in analitičnih tehnik, temelječih na inženirskih znanjih, 2. razvoj novih naprav, algoritmov, procesov in sistemov, ki podpirajo biologijo in medicino ter izboljšujejo medicinsko prakso ter skrb za zdravje« [The Whitaker Foundation, 2004]. Vsa področja so izrazito interdisciplinarna ter poleg znanj iz področja biologije, medicine, farmacije in veterine vključujejo tudi znanje s področja elektrotehnike. Uporabna so predvsem naslednja elektrotehniška znanja:
•    elektronika, kamor spada predvsem razvijanje različnih elektronskih naprav, potrebnih na tem področju; naprave so med seboj zelo različne, zahteve zanje so izrazito specifične, nekatere se proizvajajo serijsko, veliko pa je maloserijskih ali unikatnih;
•    meritve, tako elektromagnetnih in električnih signalov kakor tudi neelketričnih: tlaka, oksigenacije, stopnje glukoze … v bioloških organizmih za potrebe laboratorijskih meritev ali za implementacijo tehnologije in znanja meritev v različne naprave, npr. srčne spodbujevalnike;
1 Uvod
2
•    regulacije, saj se v bioloških organizmih srečamo s številnimi kompleksnimi regulacijskimi zankami, prav tako je znanje regulacij potrebno pri načrtovanju različnih medicinskih naprav;
•    elektromagnetika, kjer je mišljeno predvsem znanje o širjenju elektromagnetnega valovanja; zajema različne probleme, tako recimo: prevajanje električnih signalov v živčnem sistemu človeka, določanje izpostavljenosti elektromagnetnemu sevanju;
•    modeliranje, ki kot pomoč služi pri prepoznavanju in razumevanju različnih procesov, kot pripomoček ali korak pri načrtovanju naprav, posameznih sklopov in regulacijskih zank.
Cilj modeliranja je določiti obnašanje sistema. Pri tem realen sistem nadomestimo z matematičnim modelom, definiramo vhodne in izhodne spremenljivke ter izračunamo izhodne spremenljivke v odvisnosti od vhodnih spremenljivk. Matematični model, ki je v večini primerov diferencialna enačba, v primeru preprostih sistemov natančno opiše sistem, v primeru kompleksnih sistemov, kar je v biomedicini običajno, pa je matematični model poenostavljen zapis sistema. Kljub poenostavitvam dobljeni matematični model ponavadi ni analitično rešljiv, zato se je veja modeliranja izrazito razvila v zadnjem času, ko so postali dostopni računsko dovolj zmogljivi računalniki. Razvilo se je novo matematično področje, področje numeričnih metod, ki je omogočilo numerično reševanje sistema diferencialnih enačb.
Področja, kjer se modeliranje in izračun z numeričnimi metodami v biomedicini uporablja, so številna. Mednje sodijo:
•    določanje stopnje specifične absorpcije (SAR);
•    farmakokinetika in farmakodinamika;
•    računanje obremenitev v kosteh, sklepih, mišicah;
•    modeliranje delovanja organov, sistemov ali regulacijskih zank v organizmih:
o   kardiovaskularni sistem,
o   električno delovanje srca,
o   dihanje,
o   možgani,
o   mišice,
o   uravnavanje telesne temperature,
1 Uvod
3
o živčni sistem;
•    računanj e vsiljene transmembranske napetosti na celicah;
•    računanj e gostote električnega toka in porazdelitve električne poljske jakosti v tkivih;
•    populacijiski modeli;
•    modeliranje delovanja vsadljivih naprav;
•    načrtovanje terapije (npr. obsevanja v onkologiji).
1.2  Osnove modeliranja in izračuna z numeričnimi metodami
Cilj modeliranja je čim bolj natančno izračunati vrednosti spremenljivk v sistemu. Poznamo različne postopke modeliranja in izračuna z numeričnimi metodami, med najbolj razširjenimi metodami sta metoda končnih elementov in metoda končnih diferenc. Bistvo obeh metod je razdelitev celotne geometrije (območje računanja) na manjše dele, elemente, ki so v primeru končnih diferenc na premici daljice, v ravnini pravokotniki ter v prostoru kvadri, v primeru metode končnih elementov daljice na premici, v ravnini običajno trikotniki, lahko poljubni mnogokotniki ali celo liki z ukrivljenimi stranicami ter v prostoru običajno tetraedri (tristranične piramide), lahko pa poljubna geometrijska telesa, podobno kot v ravnini tudi z ukrivljenimi ploskvami. Takšni delitvi območja pravimo mreža. Vsakemu elementu v mreži nato določimo fizikalne lastnosti snovi. Na podlagi diferencialne enačbe sledi pri metodi končnih diferenc zapis enačb iskane spremenljivke za vsako vozlišče. Nastali sistem enačb pretvorimo v matrični zapis:
M-V = P,                                                                                                       (1.1)
kjer je M matrika koeficientov enačb, V vektor vrednosti iskane spremenljivke, P pa vektor robnih pogojev in virov v modelu. Dobljeni sistem enačb rešimo z eno od numeričnih metod, npr. Gaussovo eliminacijsko metodo in izračunamo vrednost iskane spremenljivke v vozliščih [Sinigoj, 1996]. Pri metodi končnih elementov je pristop k zapisu končne matrične enačbe drugačen. Kot izhodišče služi iskani spremenljivki u prirejen funkcional, iščemo pa njegov minimum. Končna matrična enačba ima enako obliko kot enačba za metodo končnih diferenc (1.1), razlikujejo se le koeficienti v matriki M in vektor robnih pogojev P.
1 Uvod
4
1.2.1 Posebnosti modeliranja v biomedicini
Modeliranje in izračun z numeričnimi metodami v biomedicini je v veliko primerih zahtevno opravilo. Kakor je napisano v prejšnjem poglavju, je potrebno za vsak element določiti lastnosti snovi v njem. Modeliranje in izračun z numeričnimi metodami uporabljamo največ na področju tehnike, kjer so lastnosti materialov znane in večinoma linearne ter izotropne, medtem ko so lastnosti tkiv v bioloških organizmih velikokrat nelinearne, anizotropne in v primeru elektromagnetnih lastnosti izrazito frekvenčno odvisne (Pavšelj, 2002). Geometrija je v tehniki večinoma sestavljena iz različnih osnovnih likov in teles, medtem ko so oblike v primeru bioloških modelov nepravilne, posledično je geometrija nesimetrična. Poleg tega je v bioloških modelih velikokrat zelo veliko razmerje med najmanjšo ter največjo dimenzijo v geometriji. Kakor bomo videli v naslednjem poglavju, predstavlja ravno veliko geometrijsko razmerje pri modeliranju zahteven problem.
1.3 Problem velikih geometrijskih razmerij
Modeliranje in izračun z numeričnimi metodami temelji na razdelitvi geometrije v elemente. Takšno delitev imenujemo mreža. V programskih paketih, namenjenih takšnemu modeliranju, se za gradnjo mreže uporabljajo posebni algoritmi za gradnjo mreže, ki morajo zagotoviti, da mreža čim bolje sledi tako robovom modela kakor tudi mejam med različnimi snovmi v modelu ter istočasno zagotoviti, da so elementi v mreži čim manj deformirani. Na primer, če mrežo sestavljajo tetraedri, bi bila mreža idealno sestavljena samo iz enakostraničnih tetraedrov, ker pa na takšen način ni mogoče slediti robovom in mejam modela, algoritem ustrezno deformira elemente, da lahko sledi meji. Od stopnje deformiranosti elementov je odvisna kvaliteta izračuna; bolj kot je element deformiran, večja je napaka izračuna. Dodatno težavo pri gradnji mreže predstavlja veliko geometrijsko razmerje (razmerje med najmanjšo in največjo dimenzijo v geometriji). Če imamo veliko geometrijo, v kateri so tudi majhna telesa, algoritem prične z gradnjo mreže v tem delu. Z oddaljevanjem mreže od začetnega mesta se velikost elementov veča. Velikokrat se dogaja, da algoritem ne more slediti takšnemu naraščanju velikosti elementov, kar privede do prekinitve gradnje mreže. Primeri, kjer so takšna velika geometrijska razmerja v biomedicini, so številni:
•    celica in celična membrana; kjer znaša debelina membrane 5 nm, premer celice pa 5 do 10 µm, razmerje je 1:1000 ali več;
1 Uvod
•    celična  suspenzija  in  celica;  kjer je razmerje odvisno predvsem od  območja suspenzije, ki jo želimo opazovati;
•    človeško telo in posamezen manjši organ; človeško telo je dolgo skoraj 200 cm, velikost posameznega organa je lahko le 1 cm ali manj, razmerje je 1:200 ali več;
•    človeško telo in elektrode; elektrode imajo premer običajno manjši kot 1 mm;
•    človeško telo in vsadki; ponavadi celotni vsadki niso večji od 5 cm, veliko pa jih ima tudi majhne elektrode: srčni spodbujevalniki, defibrilatorji, kohlearni implant …
•    človeško telo in osteosintetski material.
1.4 Možnosti za reševanje problema velikih geometrijskih razmerij
Z modeli, kjer je prisoten problem velikih geometrijskih razmerij, smo se v okviru raziskovalnega dela v Laboratoriju za biokibernetiko že srečevali. Uporabili smo različne pristope, tako smo denimo za celično suspenzijo njeno prevodnost v vektorski obliki določili na osnovi nadomestne prevodnosti celice [Pavlin in Miklavčič, 2003]. Drugi način, kako reševati problem velikih geometrijskih razmerij v modelih, je bila uporaba simetrij in antisimetrij v modelu [Susil et al., 1998]. Zaradi simetričnosti in anstisimetričnosti modela smo v primerih, ki so to omogočali, namesto izračuna v celotnem modelu izračunali le polovico, četrtino ali celo osmino osnovnega modela in nato rezultat preslikali še na ostali del geometrije. Tretji način, ki smo ga uporabljali, je bil omejitev področja računanja. V primeru lokalnega vira smo s postavitvijo izolativne meje (homogen Neumannov robni pogoj) okrog tega vira zmanjšali velik model na manjši brez prevelike napake [Šel, 2003].
Poleg številnih prednost imajo omenjeni trije načini tudi določene slabosti. Prvi, ki zahteva določanje nadomestnih prevodnosti, velja za celične suspenzije in tkivo, v primeru celične suspenzije je potrebno upoštevati še gostoto celične suspenzije. Prav tako ne omogoča opazovanja pojavov na nivoju same celice. Drugi način, kjer izkoristimo simetrije in antisimetrije modela, je mogoč samo v primeru geometrij s simetrijami. Tretji način žal ne omogoča opazovanja dogajanja večje geometrije v celoti.
1 Uvod
6
1.5 Cilji naloge
Na podlagi dosedanjih izkušenj, predstavljenih problemov teh spoznanj bomo v okviru naloge preučili možnost uporabe dveh načinov, s katerima bi lahko rešili problem velikih geometrijskih razmerij pri modeliranju in izračunu z numeričnimi metodami in jih do sedaj v Laboratoriju za biokibernetiko še nismo uporabljali.
Prva metoda je skaliranje elektrod. V Laboratoriju za biokibernetiko velikokrat rešujemo probleme, kjer kot izvor električnega polja uporabljamo elektrode, nameščene v tkivu oz. na tkivu. Elektrode so ponavadi zelo majhne (premer manj kot 1 mm [Šel et al., 2003], zato je gradnja mreže v večjem modelu zelo zahtevna. V nalogi bomo z modelom ocenili, za koliko lahko v modelu povečamo elektrode ter temu ustrezno znižamo napetost na njih, da bo rezultat še dovolj dober.
Drugi način je uporaba nelokalnega sklapljanja (non-local coupling). Pri metodi končnih elementov ali končnih diferenc na podlagi zgrajene mreže sestavimo matrično enačbo (1.1), ki jo je potrebno rešiti. Namesto da bi imeli samo eno geometrijo, bomo iz celotne geometrije izrezali tisti del, ki predstavlja največ težav za gradnjo mreže oziroma bi v njem želeli rešitev z večjo prostorsko ločljivostjo drugo, in ta del vstavili v drugo, manjšo, geometrijo. Robne pogoje na mestu rezanja bomo določili tako, da se bodo preslikali iz večje geometrije na manjšo in obratno. V vsaki geometriji bomo zgradili svojo mrežo, ki bo ločena od mreže v drugi geometriji, prav tako bodo lahko parametri pri gradnji mreže različni, matrična enačba pa bo zaradi preslikovanja robnih pogojev ena sama.
Predstavljena načina reševanja problema velikih geometrijskih razmerij pri modeliranju bioloških tkiv z metodo končnih elementov bomo v okviru naloge preizkusili na različnih modelih. Skaliranje elektrod bomo preverili na modelu dveh igelnih elektrod v homogenem ter nehomogenem tkivu, nelokalno sklapljanje pa na modelu možganov s tumorjem in vzporednimi osemigelnimi elektrodami ter na modelu celice z membrano.
2 Teorija
7
2 Teorija
2.1 Elektromagnetno polje
Splošne Maxwellove enačbe: razširjen Amperov zakon, Faradayev zakon indukcije, Gaussov stavek za magnetno polje ter Gaussov stavek za električno polje se glasijo:
riE
VxB-jU0s0 — = <u0J.                                                                                     (2.1)
dt
VxE + t = 0                                                                                                 (2.2)
V-B = 0                                                                                                          (2.3)
W-E = ^                                                                                                         (2.4)
B je gostota magnetnega pretoka, ß0 in s0 permeabilnost in dielektričnost praznega prostora,
E električna poljska jakost, J volumska gostota električnega toka ter p prostorska gostota naboja. V elektrotehniki se pogosteje uporabljajo Maxwellove enačbe, ki upoštevajo polarizacijo in magnetizacijo snovi. Tako dobimo dve novi spremenljivki, gostoto električnega pretoka D in magnetno poljsko jakost H :
H = B   ,                                                                                                       (2.5)
D = LrL0E                                                                                                        (2.6)
kjer sta /ur in sr relativna permeabilnost in relativna dielektričnost snovi. Po tako vpeljanih spremenljivkah se sistem Maxwellovih enačb glasi:
2 Teorija
8
VxH-— = Jprosti.                                                                                         (2.7) dt
VxE += 0                                                                                                 (2.8) dt
V-B = 0                                                                                                          (2.9)
y.D = pprosti                                                                                                                   (2.10)
Z  Jprosti   in  Pprosti   sta označena volumska gostota električnega toka prostih nabojev in
volumska gostota prostih nabojev, to je tistih nabojev v snovi, ki niso vezani in se lahko premikajo. Takšen sistem enačb je v elektrotehniki bolj pogost od prvega, predstavljenega v enačbah (2.1) do (2.4).
Poleg Maxwellovih enačb so za poznavanje elektromagnetnega polja pomembne še naslednje enačbe:
V-J + ^ = 0,                                                                                               (2.11)
dt
ki ga poznamo kot zakon o ohranitvi nabojev oz. kontinuitetni zakon za električni tok,
E = -WV,                                                                                                      (2.12)
ki govori o tem, da je električna poljska jakost gradient električnega potenciala V,
B = VxA,                                                                                      (2.13)
ki definira magnetni vektorski potencial A in Ohmov zakon:
J = yE.                                                                                                        (2.14)
Specifična električna prevodnost snovi rje v splošnem tenzor:
2 Teorija
9
y =
I   XX          i   .
xy
Yx
yx          yy          yz
i zx        i zy        i z
(2.15)
vendar se da s pravilno orientacijo za večino snovi zapisati tako, da so vsi izven diagonalni elementi enaki 0. Sicer pa je veliko snovi izotropnih, kar pomeni, da je njihova specifična prevodnost v vseh smereh enaka in jo torej lahko zapišemo s skalarjem.
Pri prehajanju elektromagnetnega polja preko meje dveh snovi je zelo pomembno poznavanje pravil prehajanja omenjenega polja preko te meje. Z indeksi 1 označimo veličine v prvi snovi in z indeksi 2 veličine v drugi snovi, oznaka t v indeksu pomeni tangencialno komponento veličine glede na ravnino meje, oznaka n pa normalno komponento veličine glede na ravnino meje. Iz prvih dveh Maxwellovih enačb izhajata prestopna pogoja za poljski jakosti:
H t 2 -H t1
K
prosti
(2.16)
Kprosti je površinska gostota električnega toka na meji.
F   -F   -0
(2.17)
Iz drugih dveh Maxwellovih enačb izhajata prestopna pogoja za gostoti pretoka:
Bn2-Bn1=0,
(2.18)
D   -D 1 =<7       .
n2           m           prosti
(2.19)
?prosti je površinska gostota naboja na meji.
2 Teorija
10
2.1.1 Tokovno polje
Določanje vseh veličin elektromagnetnega polja je računsko zahtevno in zamudno. Programski paketi za modeliranje in numerično reševanje to omogočajo le deloma, večinoma le za harmonska vzbujanja v stacionarnih razmerah. Velikokrat lahko računanje poenostavimo tako, da namesto vseh veličin določimo le tiste, ki so pomembne. Tako je na primer mogoče stimulacijo bioloških tkiv z enosmernim električnim tokom opisati s tokovnim poljem. Tokovno polje je vezano na prevodne snovi, kjer velja Ohmov zakon. Od sistema štirih Maxwellovih enačb ostaneta le dve poenostavljeni enačbi:
V-E = ^                                                                                                       (2.20)
VxE = 0                                                                                                       (2.21)
Poleg tega veljata tudi (2.14) Ohmov zakon in (2.11) zakon o ohranitvi nabojev oz. kontinuitetni zakon za naboj. Ko v enačbo (2.11) vstavimo enačbi (2.14) ter (2.12), dobimo enačbo:
v-(r-(-w)) = -.                                                                (2.22)
Če velja, da je specifična električna prevodnost snovi  y izotropna in homogena (y je konstanta), dobimo Poissonovo diferencialno enačbo:
AV = ^.                                                                                                      (2.23)
ydt
V stacionarnih razmerah, ko so vrednosti časovnih odvodov enake 0, se enačba poenostavi v homogeno Poissonovo oz. Laplaceovo diferencialno enačbo:
AV = 0.
(2.24)
2 Teorija
11
2.2 Modeliranje in izračun z numeričnimi metodami
Uporaba modeliranja se veča. Vzrok temu je predvsem v razvoju računskih zmogljivosti računalnikov ter njihovi dostopnosti, po drugi strani pa tudi v večji ponudbi programskih paketov za modeliranje in izračun z numeričnimi metodami ter vse večja uporabnost le-teh.
Pri modeliranju sistem nadomestimo z matematičnim modelom, ki predstavlja sistem v obliki diferencialnih enačb. Kljub poenostavitvam dobljene diferencialne enačbe ponavadi niso analitično rešljive. Pred uporabo digitalnih računalnikov so se za reševanje takšnega sistema diferencialnih enačb uporabljali analogni računalniki. Analogni računalniki so vsebovali različne osnovne gradnike, kot na primer različni viri, množilnik, seštevalnik, diferenciator, integrator …, ki jih je bilo potrebno med seboj pravilno povezati glede na diferencialne enačbe [Zupančič et. al., 1995]. Danes so analogne računalnike nadomestili digitalni računalniki, ki so postali dovolj zmogljivi, da omogočajo numerično reševanje sistema diferencialnih enačb.
Modeliranje in izračun z numeričnimi metodami prinaša številne prednosti, saj lahko na primer preizkušamo občutljivosti sistemov na posamezne parametre, kar je na samih sistemih ponavadi veliko dražje, mnogokrat pa to zaradi različnih vzrokov nemogoče (varnost, čas, zdravje …); pri načrtovanju sistema se lahko opiramo na podatke, pridobljene iz odzivov modela, in podobno.
V okviru Laboratorija z biokibernetiko uporabljamo tri programske pakete za modeliranje in izračun z numeričnimi metodami: FEMLAB (Comsol, Švedska), MAXWELL (Ansoft, ZDA) in EMAS (Ansoft, ZDA). Vsi trije programski paketi temeljijo na metodi končnih elementov, vendar se za uporabnika med seboj precej razlikujejo.
2 Teorija
12
2.2.1 Metoda končnih diferenc
Za razumevanje numeričnih metod je najprimernejša metoda končnih diferenc. Za tokovno polje temelji zapis metode na enačbi (2.22)
V-(y(-VV)) = -^
dt
ki se v stacionarnih razmerah poenostavi v enačbo (2.25):
V.(r(-VV)) = 0.                                                                                         (2.25)
Za nek majhen volumen W (oz. za zapis v integralski obliki) bo veljalo:
jW(-f7V)dw = 0.                                                                                        (2.26)
W
Če uporabimo Gaussov izrek in nadomestimo volumski integral divergence s ploskovnim integralom po sklenjeni ploskvi, dobimo iz enačbe (2.26):
-j)f7Vda = 0                                                                                               (2.27)
A
ali drugače zapisano:
Ay^V da = 0.                                                                                             (2.28)
A        dn
Enačba (2.28) je izhodiščna enačba za računanje porazdelitve električnega potenciala v primeru tokovnega polja z metodo končnih diferenc.
2 Teorija
13
C
1 BC	hC	I AC
hB    Ä		hA
^		
1 BD	hD	1 AD
D
Slika 3.1: Vozlišče v kvadratni mreži v ravnini. Z A, B, C in D so označena sosednja vozlišča v mreži, h označuje razdalje do sosednjih vozlišč, xpa specifične prevodnosti.
Pri izračunu z metodo končnih diferenc je potrebno želeno območje prekriti z daljicami na premici, s pravokotniki v ravnini in kvadri v prostoru. V ravnini nadomestimo integral v enačbi (2.28) za eno vozlišče, predstavljeno na sliki 3.1, z vsoto posameznih prispevkov:
&vdV           (v    hCd + v    hDdIVA   V 0    (v    hCd         hDdVB-V0
Ar-da^ — ^—) — -^ — ^—) — -^ hB
_fv   hBd+v  hAdVC~V 0  L   hBd | ,,   hAdIM ^ 2 +rAC 2 j   hC     ^ 2 +rAD 2 j   hD
kjer je d debelina plasti. Z A, B, C in D so označena sosednja vozlišča, h pa označuje razdalje do sosednjih vozlišč. Takšno enačbo zapišemo v vseh vozliščih ter iz njih sestavimo matrični zapis (1.1):
M -V = P ,
kjer je M matrika koeficientov enačb, V vektor vrednosti iskane spremenljivke, P pa vektor robnih pogojev in virov v modelu. Dobljen sistem enačb rešimo z eno od numeričnih metod in izračunamo vrednost iskane spremenljivke v vozliščih. Podobno enačbo kot (2.29) dobimo tudi v primeru prostora, le da imamo na desni strani enačbe namesto sedanjih štirih šest členov.
2 Teorija
14
2.2.2 Metoda končnih elementov
Kakor smo že prej napisali, je končna matrika za reševanje pri metodi končnih elementov podobna kot pri metodi končnih diferenc, prav tako je osnova obeh metod razdelitev prostora v manjša območja, le da mreža pri metodi končnih diferenc temelji na kvadratnih oblikah, pri metodi končnih elementov pa običajno temelji delitev na trikotnih oblikah, lahko pa na poljubnih mnogokotnikih, celo takšnih z ukrivljenimi robovi ali stranicami. Vendar se metodi razlikujeta po pristopu k zapisu enačb, saj pri metodi končnih diferenc zapis temelji na Poissonovi enačbi in robnih pogojih, pri metodi končnih elementov pa zapis temelji na minimizaciji potencialu prirejenega funkcionala.
Prvi korak pri metodi končnih elementov je razdelitev območja računanja v v elemente, kar imenujemo tudi gradnja mreže. Na premici predstavljajo posamezne elemente daljice, ki so med seboj lahko različno dolge, ne smejo se prekrivati, hkrati pa morajo pokriti celotno območje računanja v. Podobno je v ravnini, ki jo ponavadi prekrijemo s trikotniki in prostoru, kjer za elemente ponavadi vzamemo tetraedre. Pomembno je, da elementi s čim manjšim odstopanjem sledijo robovom modela in mejam med različnimi snovmi v modelu ter so obenem tako majhni, da lahko v posameznem elementu lastnosti snovi opišemo kot konstantne. Sledi določanje normiranih lokalnih koordinat v vsakem elementu. V primeru tetraedra, katerega oglišča označimo z A, B, C in D, uporabimo piramidne koordinate gA, ÇB, ÇeC in Ž d, ki imajo izhodišča v ogliščih in naraščajo v smeri pravokotnic do nasprotnih trikotnih stranic tetraedra do vrednosti 1. Določanju elementov sledi določanje vozlov. Najpreprostejši je primer, ko so vozli postavljeni le v krajišča daljic oziroma v oglišča trikotnikov ali tetraedrov, lahko pa so tudi na sredini daljic oziroma stranic ali robov. Prvi primer imenujemo linearni elementi, drugega kvadratni, če pa bi uporabili še dodatne vmesne vozle, bi naraščala tudi potenca elementov.
Imejmo na območju računanja v dano funkcijo y/ = yix, y, z), ki na ogliščih elementa zavzame vrednosti y/A, Vb, Vc in yD. Aproksimacijo funkcije y/ v enem elementu, y/e definiramo tako, da v ogliščih zavzame enake vrednosti kot funkcija y/, kar je najenostavneje doseči z linearno funkcijo \j/e = otx + ßy + /z + &, ki jo izrazimo kot funkcijo normiranih koordinat elementa e:
2 Teorija
15
fe(x,y,z) = <pe(šeA,šeB,še,šeD) = <pe.                                                                  (2.30)
Od tod tudi ime linearni elementi. Ob uporabi dodatnih vmesnih vozlišč bi namesto linearne funkcije za zapis \j/e uporabili kvadratno in elemente imenovali kvadratne. Aproksimacijo cpe sestavimo iz vsote štirih linearnih funkcij:
D
če = Lm1-L).                                                                      (2.31)
T=A
Taka aproksimativna funkcija cpe v vseh vozliščih zavzame enake vrednosti kot funkcija y/, kar tudi želimo. Za zapis aproskimacije funkcije y/v vseh elementih modela, if/ oziroma q> definiramo funkcijo
e      öe za točke v elementu
ç»l=\                                                                                                            (2.32)
0 drugje
ter iz vsote takšnih aproksimacijskih funkcij določimo aproksimacijo celotnega območja računanja
<P = E>0,                                                                                                     (2.33)
kjer je E število vseh elementov. Z indeksom i označimo vse vozle v mreži: i = 1, 2, …, n ter s \j/i vrednosti funkcije y/v vozlih in zapišimo q> kot dvojno vsoto:
En                   n          E               n
v=Ž Z viuie =L ^Z uie = Z ^iui ,                                                           (2.34)
e=1  i=1                      i=1          e=1                 i=1
kjer je
2 Teorija
16
e     \1-<Ee velementu u =\      hi                  .                                                                                  (2.35)
[0        drugje
Zadnji del enačbe, vsota uie po vseh elementih, se v enem vozlu reducira le na tiste elemente, v katerih leži to vozlišče. To vsoto imenujmo ui, to je vsota vseh tistih funkcij {1-f i}, ki linearno upadajo iz i-tega vozlišča proti sosednjim vozliščem.
Ko imamo že določeno mrežo elementov v v, sledi določanje iskanih veličin, v primeru tokovnega ali elektrostatičnega polja je to električni potencial V. Nastavek za potencial zapišemo v obliki
V = LCiui,                                                                                                   (2.36)
i
torej kot linearno kombinacijo funkcij {ui} s koeficienti {Ci}, ki jih ne poznamo. Koeficienti Ci zavzamejo v vozliščih aproksimativne vrednosti električnega potenciala, Ci = Vi iz česar sledi:
i
Razširimo zadnjo enačbo z vpeljavo Dirichletovih robnih pogojev za tista vozlišča, kjer takšni robni pogoji veljajo, denimo da je m takšnih vozlišč.
mn
1 V 0ku0k+   1
k=                        i=
iti                                            it
V = YV 0ku0k +LV iui                                                                               (2.38)
Funkcije {u0k} so definirane enako kot funkcije {ui}, {V0k} pa so znane vrednosti potencialov v teh vozlih. Z uporabo Rayleigh-Ritzovega postopka [Sinigoj, 1996] za minimizacijo funkcionala iz enačbe (2.38) dobimo za elektrostatično polje enačbo:
n   f                             \             f                        m                                   \
2  jXvuiVuj) \V j=\\ Pprostiui-e^jV^Vui-Vu0 k) dv+ j suih 'da,         (2.39)
j=1 V v'                            J           v\                        k=1v'                                  J          An'
kjer je v' območje računanja, ki ga opisujejo elementi in odstopa od v zaradi oglatosti robov, AN' je meja v', s je relativna dielektričnost snovi, pprosti prostorska gostota naboja v snovi ter h '
2 Teorija
17
odsekoma konstantna aproksimacija normalnega odvoda. Enačba (2.39) za tokovno polje je podobna, le da relativno dielektričnost snovi s v enačbah nadomesti specifična prevodnost snovi fter prostorsko gostoto naboja pprosti nadomesti odvod prostorske gostote naboja po času, ki pa je za stacionarne razmere enak 0:
n   f                             \             f        m                                   \
S  JV(Vui • Vuj) \V j = j] -rL jV0k (Vui • Vu0k) dv+ j yuih 'da,                     (2.40)
j=1 V                         J          v\      k=1v'                              J         AN'
Po reševanju vseh integralov v enačbi (2.40) [Sinigoj,  1996] dobimo sistem enačb za aproksimativni potencial v vsakem elementu Vi v naslednji obliki:
Y1 mijVj=pi,   i = 1,2,...,n,                                                                            (2.41)
kjer so
ter
I             m                                                     \
pi=\\ -rH\V0 k{Vui-Vu0k) dv+ j yuih'da
J               AN'
V matrični obliki bi enačbo (2.41) zapisali kot
M -V = P
kar pa je pravzaprav enačba (1.1). Rešitev sistema enačb (2.41) vstavimo v enačbo (2.38) ter tako določimo iskan aproksimativni potencial V .
v'
3 Reševanje problema velikih geometrijskih razmerij
18
3 Reševanje problema velikih geometrijskih razmerij
3.1   Uporaba skaliranih elektrod
V okviru različnih izračunov se v Laboratoriju za biokibernetiko velikokrat srečamo s problemom, ko želimo določiti električno polje v tkivu, kjer kot vir električnega polja uporabljamo igelne elektrode. Igelne elektrode so ponavadi debele manj kot 1 mm, npr. igle, uporabljane pri Cliniporatorju, imajo premer 2r = 0.6 oz. 0.7 mm. Velikost takšnega modela omejuje algoritem za gradnjo mreže v različnih komercialnih programskih paketih za računanje z metodo končnih elementov, zato je nemogoče določiti električno polje v večji oddaljenosti od elektrod. Če namreč želimo izračun v bližini elektrod, lahko brez prevelike napake v primerni oddaljenosti od elektrod omejimo naš prostor z homogenim Neumannovim robnim pogojem, kar pomeni, da tam naš model izoliramo. Primerno oddaljenost za določitev meje modela najlažje določimo s preizkušanjem, splošno pravilo pa je, naj bo nekajkrat večja od največje dimenzije objekta, v katerem nas porazdelitev električne poljske jakosti ali električnega potenciala zanima. S tem pristopom pa ne moremo izračunati porazdelitve v večji oddaljenosti od elektrod. Možna rešitev bi bila uporaba nadomestnih skaliranih elektrod, ki bi imele premer rs za faktor skaliranja k večji:
k = 100-( rs—1).                                                                                               (3.1)
r
Ker bi bile nadomestne elektrode debelejše, bi bila gradnja mreže enostavnejša, seveda pa bi bilo potrebno temu ustrezno tudi znižati priključeno napetost na njih. Za nehomogen model na sliki 3.1 smo preverili, kolikšno je lahko skaliranje elektrod, ne da bi bila pri tem napaka pri izračunu prevelika. Kot izhodišče smo uporabili dve igelni elektrodi, premera 0.7 mm in dolžine 40 mm. Elektrodi imata električni potencial V = ±100 V. Igli sta bili postavljeni v središče nehomogene kocke z robom 200 mm, sestavljene iz dveh polovic, in sicer prve s
3 Reševanje problema velikih geometrijskih razmerij
19
prevodnostjo ?1 = 0.1 Sm-1 ter druge s prevodnostjo ?2 = 1 Sm-1. Razmerje prevodnosti v tkivih je običajno manjše od zgoraj omenjenega razmerja 10, saj velja za: mišico transferzalno ? = 0.04-0.14 Sm-1, mišico longitudinalno ? = 0.3-0.8 Sm-1, nižje plasti kože (epidermis, dermis) ? = 0.2 Sm-1, nižjo prevodnost imajo le maščoba ? = 0.02-0.04 Sm-1 ter kosti ? = 0.006 Sm-1 [Pavšelj, 2002]. Za reševanje sistema enačb v FEMLAB-u smo uporabili direktni linearni algoritem ter izračune opravili za naslednje razdalje med elektrodama d = 4, 6, 8 in 10 mm, oziroma razdalje med osema elektrod d0 = 4.7, 6.7, 8.7 ter 10.7 mm. Opazovali smo ekvipotencialne ploskve pri električnem potencialu Vs = ±25 V, ±50 V ter ±75 V, kar ustreza 1, 1 ter 3 električnega potenciala elektrod V. Nadomestne elektrode morajo s čim manjšim odstopanjem nadomestiti ekvipotencialne ploskve, tako da je popačenje električnega polja v oddaljenosti majhno. Na podlagi volumna, ki ga objamejo ekvipotencialne ploskve, smo določili, kolikšen je premer nadomestnih elektrod ter iz tega izračunali faktor skaliranja premera oziroma polmera k za posamezno razmerje med električnim potencialom elektrod ter električnim potencialom ekvipotencialnih ploskev.
Slika 3.1: Nehomogen prostor dimezije 200 × 200 × 200 mm z dvema snovema, prva ima prevodnost ?1 = 0.1 Sm-1, druga pa ?2 = 1 Sm-1. Razdalja med elektrodama d na sliki je 4 mm, premer 0.7 mm ter dolžina 40 mm.
3 Reševanje problema velikih geometrijskih razmerij
20
Za ugotavljanje napake zaradi uporabe skaliranih elektrod smo rezultate izhodiščnega modela primerjali z rezultati modela s skaliranimi elektrodami ter določili odstopanje v porazdelitvi električne poljske jakosti in električnega potenciala.
3.2 Uporaba nelokalnega sklapljanja v FEMLAB-u
Gradnja modela v programskem paketu za modeliranje z metodo končnih elementov FEMLAB temelji na korakih, opisanih v poglavju 2.2.2, z nekaj posebnostmi:
•    določitev oz. risanje geometrije;
•    izgradnja mreže;
•    določanje fizikalne narave problema (application mode), kajti parametri v naslednjih korakih so odvisni od fizikalne narave problema;
•    določanje konstant, spremenljivk in sklapljanja;
•    določanje robnih pogojev;
•    določanje lastnosti snovi;
•    izbor algoritmov in parametrov za reševanje ter reševanje samo in
•    predstavitev rezultatov.
V FEMLAB-u smo zgradili dva modela, ki sta za gradnjo mreže zelo zahtevna. Prvi model je model človeških možganov s tumorjem. Povzet je po doktorski nalogi Davorke Šel [Šel, 2003]. Glavna težava reševanja tega problema je bila gradnja mreže v modelu po vstavitvi vzporednih osemigelnih elektrod. Drug model je model celice z membrano. Zaradi tankosti membrane v primerjavi s premerom celice je nemogoče zgraditi mrežo v modelu celotne celice. Oba problema smo reševali z uporabo nelokalnega sklapljanja.
Pri nelokalnem sklapljanju smo iz celotne geometrije izrezali tisti predel, ki je predstavljal težavo pri gradnji mreže oziroma smo želeli v njem prostorsko bolj podroben ter natančen rezultat. Izrezani del smo vstavili v drugo geometrijo. Tako smo dobili dve geometriji v enem modelu: celotno brez izrezanega dela, imenovali smo jo velika, ter izrezani del, ki smo ga imenovali mala geometrija. FEMLAB namreč omogoča, da lahko v okviru enega modela določimo več geometrij.
3 Reševanje problema velikih geometrijskih razmerij
21
Kot fizikalno naravo modela smo izbrali stacionarno tokovno polje. Tako v mali kot v veliki geometriji je bilo z robnimi pogoji na meji rezanja, ki je postala rob, potrebno zagotoviti, da se bo meja rezanja obnašala kakor meja dveh prevodnih medijev. Na meji dveh prevodnih medijev velja zaradi zakona o ohranitvi nabojev za stacionarno tokovno polje (2.11) naslednji prestopni pogoj za gostoto električnega toka:
Jn2-Jn1=0.                                                                                                   (3.2)
Zaradi nevrtinčnosti električne poljske jakosti (2.2) velja za električno poljsko jakost naslednji prestopni pogoj:
Et2-Et1=0.                                                                                                  (2.17)
Poleg izpolnjevanja mejnih pogojev je bilo potrebno na meji rezanja zagotoviti še zveznost električnega potenciala V, kar je v primeru mejnih pogojev samoumevno, v primeru robnih pogojev pa ne. Takšni robni pogoji so torej sestavljeni robni pogoji, saj moramo na robu poznati električni potencial V ter njegov odvod, ki je sestavljen iz dveh komponent:
Et2=Et1.                                                                                                        (3.3)
in
En2=LE n1.                                                                                                     (3.4)
I2
FEMLAB omogoča reševanje takšnega problema, saj omogoča definiranje več vrst spremenljivk sklapljanja (coupling variables) [Comsol AB., 2001]. To so spremenljivke, ki so lahko definirane v eni geometriji in jih uporabimo tako v tej kot tudi v drugi geometriji znotraj istega modela. Omogočajo nam preslikovanje med različnimi dimenzijami v geometriji, na primer integral določene spremenljivke po volumnu v eni geometriji služi kot robni pogoj v drugi geometriji. Z eno od takšnih spremenljivk sklapljanja (extrusion coupling variable) smo na meji rezanja določili sklopljen robni pogoj. Pri tem ni potrebno določiti sestavljenega robnega pogoja, ampak je potrebno na meji rezanja določiti le spremenljivko sklapljanja Vs, je definirana v večji geometriji in cilja na manjšo geometrijo. Za vrednost
3 Reševanje problema velikih geometrijskih razmerij
22
spremenljivke sklapljanja zahtevamo, da je Vs = V, spremenljivka V označuje električni potencial v večji geometriji. Nato je potrebno v manjši geometriji na meji rezanja določiti robni pogoj: V2 = Vs, na meji rezanja v večji geometriji pa homogene Neumannove robne pogoje. Spremenljivka V2 predstavlja električni potencial v manjši geometriji. Za algoritem reševanja takšno sklapljanje pomeni, da rob, na katerem je definirana spremenljivka sklapljanja, ter rob, na katerega cilja spremenljivka sklapljanja, ni več rob, ampak samo meja dveh prevodnikov. Zato algoritem zgradi eno samo matriko za reševanje, ki jo reši enako, kakor da gre za eno geometrijo. Pri tem je potrebno upoštevati še, da je matrika, ki jo je potrebno rešiti, običajno redka. Vrednost nekega vozlišča je namreč odvisna le od vrednosti sosednjih vozlišč, vsi ostali elementi v vrstici matrike pa so enaki nič. Zaradi nelokalnega sklapljanja postane ta matrika vse manj redka, zaradi česar se poveča raba pomnilnika med izračunom ter se podaljša čas reševanja [Comsol AB., 2001].
3.2.1 Model možganov s tumorjem
V programskem paketu FEMLAB je mogoče geometrijo narisati v samem paketu, mogoče pa jo je tudi uvoziti. Ena izmed možnosti je uvoz slik, iz katerih nato sestavimo geometrijo. Ta metoda je bila uporabljena pri določanju geometrije modela človeških možganov s tumorjem in ga je v okviru doktorata opravila Davorka Šel [Šel, 2003], kjer je tudi podrobneje opisan postopek obdelave slik. Za gradnjo geometrije je bilo uporabljenih 20 sivinskih CT slik velikosti 512 × 512 točk, slikanih z razmakom 5 mm. Ker bi bila geometrija preveč obsežna, je za model človeških možganov služilo le 10 slik, torej vsaka druga. Slike so bile uvožene v MATLAB (The MathWorks, ZDA), kjer so slike predstavljene kot matrike z enako velikostjo kot je velikost slik, vrednosti v matriki pa ustrezajo vrednosti barve v tej točki. Z MATLAB-ovo funkcijo »contour« so bile določene krivulje, kjer barva v sliki ustreza izbranemu sivinskemu pragu. Ker je bilo v posamezni sliki več takšnih krivulj, je bilo potrebno s posebnim algoritmom izbrati tisto, ki predstavlja zunanjo mejo možganov. Bistvo algoritma je bilo določanje velikosti sklenjene krivulje. Največja krivulja je bila krivulja na robu lobanje, druga največja pa je bila iskana krivulja na robu možganov. Vsako krivuljo je bilo potrebno nadomestiti z nepravilnim mnogokotnikom z enakim številom oglišč. V zadnjem koraku oglišča med seboj povežemo in tako nastane celotna geometrija. Podoben postopek je bil uporabljen za izgradnjo modela tumorja, le da so bile v tem primeru uporabljene vse slike, ki so vsebovale tumor, in ne le vsaka druga. V model so vključene tudi vzporedne osemigelne elektrode, kakršne se uporabljajo pri Cliniporatorju. Dolge so 40 mm, premer znaša 0.6 oz.
3 Reševanje problema velikih geometrijskih razmerij
23
0.7 mm, razdalja med vrstama elektrod je 8 mm, razdalja med elektrodama v isti vrsti pa je 3.2 mm. Model z elektrodami premera 0.7 mm je prikazan na sliki 3.2.
Slika 3.2: Model človeških možganov s tumorjem in vzporednimi osemigelnimi elektrodami.
Gradnja mreže v takšnem modelu je zahtevna, saj [Šel, 2003]:
•    geometrija vsebuje veliko ukrivljenih objektov;
•    se ukrivljeni objekti zelo razlikujejo po velikosti;
•    razdalja med posameznimi mejami je zelo majhna (tumor leži na robu možganov).
Kljub avtomatizirani gradnji mreže v FEMLAB-u je bil proces izgradnje mreže dolgotrajen, saj je bilo potrebno skrbno iskanje primernih parametrov algoritma za gradnjo mreže. Ker bi bila gradnja mreže v modelu z elektrodami prezahtevna, smo jih iz modela izrezali. Zaradi primerjave rezultatov z modelom z nelokalnim sklapljanjem, smo vse lastnosti snovi in robne pogoje določili v obeh modelih enako. Tako smo kot fizikalno naravo izbrali tokovno polje, vse robove modela smo nastavili na izolativne, razen robove, kjer bi se elektrode dotikale tkiva ali tumorja. Robovom v eni vrsti smo določilo električni potencial -300 V, robovom v drugi vrsti pa +300 V. Za prevodnost tkiva ?tkivo smo izbrali 0.1 Sm-1, za prevodnost tumorja ?tumor pa 0.3 Sm-1. Za reševanje sistema enačb v FEMLAB-u smo uporabili iterativni linearni algoritem z nepopolno L-U dekompozicijo za začetno transformacijo (preconditioning) ter Good Broyden metodo za rešitev sistema enačb.
3 Reševanje problema velikih geometrijskih razmerij
24
3.2.2 Model možganov s tumorjem in nelokalnim sklapljanjem
Zaradi zahtevnosti geometrije smo se odločili model možganov s tumorjem rešiti tudi z nelokalnim sklapljanjem. Iz geometrije možganov s tumorjem in vzporednimi osemigelnimi elektodami smo izrezali predel tumorja z elektrodami in izrezani del vstavili v drugo geometrijo v istem modelu.
Pri določanju meje rezanja je bilo potrebno le-to izbrati tako, da meja rezanja ni potekala v istih višinah, kot ležijo robovi modela, bolj natančno, da meja rezanja ni potekala v istih višinah, kot so bile osnovne konture, pridobljene iz CT slik in iz katerih je zgrajena geometrija modela. Če tega nismo upoštevali, je bila zaradi majhnih razdalj med mejo rezanja in robom modela tam mreža zelo gosta, včasih celo avtomatski algoritem za gradnjo mreže lete ni uspel zgraditi. Obe geometriji novega modela sta predstavljeni na sliki 3.3. Kot fizikalno naravo modela smo izbrali tokovno polje. Na meji rezanja smo s pomočjo spremenljivke sklapljanje določili sklopljen robni pogoj. Na vse ostale robove, razen na konce elektrod, smo nastavili izolativne robne pogoje. Konce elektrod v eni vrsti smo postavili na električni potencial -300 V, v drugi vrsti pa na +300 V. Prevodnosti so bile enake kot v prejšnjem primeru, in sicer ?tkivo = 0.1 Sm-1, ?tumor = 0.3 Sm-1 ter ?elektrod = 4×107 Sm-1. Za reševanje sistema enačb v FEMLAB-u smo uporabili direktni linearni algoritem.
3 Reševanje problema velikih geometrijskih razmerij
25
Slika 3.3: Model človeških možganov s tumorjem in vzporednimi osemigelnimi elektrodami. Na zgornji sliki je predstavljena prva, velika geometrija modela: možgani z izrezanim predelom tumorja in elektrod. Na spodnji
sliki je predstavljen izrezani del ali mala geometrija.
3 Reševanje problema velikih geometrijskih razmerij
26
3.2.3 Model celice z nelokalnim sklapljanjem
Eno izmed področji raziskovanja Laboratorija za biokibernetiko je elektroporacija. Proces temelji na uporabi električnega polja visokih jakosti, ki povzroči vsiljeno transmembransko napetost na celici. Posledica vsiljene transmembranske napetosti je lahko povečana prepustnost celične membrane. Eden izmed najpomembnejših parametrov povečanja prepustnosti membrane je velikost vsiljene transmembranske napetosti. Analitično jo je mogoče določiti za eno samo celico sferoidne oblike [Kotnik in Miklavčič, 2000, Valič et al., 2003] v neskončnem prostoru, pa še to le za neprevodno membrano. Modeliranje celice z membrano in izračun z numeričnimi metodami je nemogoč zaradi majhne debeline membrane, ki znaša nekaj nm, v primerjavi s premerom celice, ki znaša nekaj µm. Želeli smo preveriti, kakšne možnosti ponuja metoda nelokalnega sklapljanja, zato smo zgradili model celice v FEMLAB-u z nelokalnim sklapljanjem.
Slika 3.4: Geometrija modela celice. Kocka z robom 40 µm predstavlja okolico celice, krogla pa celico s polmerom 10 µm. Dve piki na skrajnih koncih krogle po osi y predstavljata dva izrezana dela, ki sta vstavljena v
novi geometriji, predstavljeni na sliki 3.5.
3 Reševanje problema velikih geometrijskih razmerij
27
Osnova geometrije je krogla polmera 10 µm v središču kocke z robom 40 µm. Krogla predstavlja celico, kocka pa okolico. Geometrija modela je prikazana na sliki 3.4. V model smo želeli vključiti tudi membrano, ki je običajno debela nekaj nm, v modelu smo za debelino membrane izbrali 5 nm. Seveda tako tanke membrane z metodo nelokalnega sklapljanja ne moremo modelirati na vsej celici. Na obeh skrajnih koncih celice po osi y, kjer so pričakovane največje vrednosti električne poljske jakosti, smo izrezali dva manjša dela celice osnovne geometrije, velika 100 nm × 30 nm × 100 nm. Na sliki 3.4 se vidita izrezana dela kot dve majhni piki. Oba izrezana dela smo vstavili v novi geometriji, imenujmo tisto, ki leži na negativni strani osi y, minus, drugo pa plus. V geometriji minus in plus smo dodali še izrezka krogle, ki ima polmer za 5 nm večji od polmera krogle, ki predstavlja celico, torej 10.005 µm. Tako smo dobili v obeh manjših geometrijah model dela celice z membrano, kakor je predstavljeno na sliki 3.5.
Slika 3.5: Geometrija modela celice. Levo je predstavljena geometrija plus (izrezani del celice na pozitivnem koncu osi y), desno je predstavljena geometrija minus (izrezani del celice na negativnem koncu osi y).
Žal algoritem za gradnjo mreže le-te v takšnem modelu ni uspel zgraditi. Problem se je pojavil v geometrijah plus in minus, čeprav sta geometriji zelo preprosti. Predvidevamo, da je razlog za to v načinu definiranja geometrije v FEMLAB-u. Kljub temu, da je lokalna geometrija enostavna, je zapisana kot funkcija osnovnih geometrij: geometrijo plus tako sestavljajo kvader z dimenzijami 100 nm × 30 nm × 100 nm (okolica), presek krogle s polmerom 10.005 nm in kvadrom (celična membrana) ter presek krogle s polmerom 10 nm in
3 Reševanje problema velikih geometrijskih razmerij
28
kvadrom (notranjost celice). Namesto takšnih geometrij plus in minus smo uporabili geometriji na sliki 3.6. Namesto krogelnih odsekov smo geometriji sestavili iz treh kvadrov: 100 nm × 30 nm × 100 nm (okolica), 100 nm × 20 nm × 100 nm (celična membrana) ter 100 nm × 15 nm × 100 nm (notranjost celice). Napaka zaradi uporabe nadomestne geometrije je zanemarljiva, kar je razvidno iz primerjave slik 3.5 ter 3.6. V takšni geometriji modela gradnja mreže ni več predstavljala težave.
Slika 3.6: Geometrija modela celice. Levo je predstavljena nova geometrija plus, desno pa nova geometrija minus. Namesto izrezov iz krogel sta sestavljeni iz kvadrov. Napaka zaradi neukrivljenosti ploskev je zelo
majhna in je na sliki ni mogoče opaziti.
Kot fizikalno naravo modela smo izbrali stacionarno tokovno polje. Na meji rezanja smo tako v veliki geometriji kot tudi v geometrijah plus in minus s pomočjo spremenljivke sklapljanja določili sklopljen robni pogoj. Na ostalih robovih velike geometrije smo določili izolativne robne pogoje, razen na stranici, vzporedni z ravnino xz, na kateri smo postavili na električni potencial -100 V ter +100 V. Posebno težavo je predstavljala meja med kroglo in kocko, ki predstavlja celično membrano. V FEMLAB-u je mogoče določati tudi robne pogoje na meji dveh snovi (notranji robni pogoj). Vendar je nemogoče nastaviti takšen robni pogoj, ki bi predstavljal prevodno plast v smeri normale na mejo. Namesto tega smo v veliki geometriji uporabili dve fizikalni naravi. FEMLAB namreč omogoča, da v eni geometriji hkrati računamo dve fizikalni naravi: tako bi lahko denimo hkrati računali prevajanje toka in prevajanje toplote, in to dvoje tudi povežemo: ohmske izgube pri prevajanju toka služijo kot vir segrevanja pri prevajanju toplote. Za vsako fizikalno naravo posebej lahko nastavimo
3 Reševanje problema velikih geometrijskih razmerij
29
robne pogoje in lastnosti snovi, prav tako pa lahko za vsako fizikalno naravo posebej nastavimo, kateri objekti v geometriji so aktivni in kateri ne. Tako smo v veliki geometriji določili dve fizikalni naravi prevajanja toka, v eni fizikalni naravi je bila aktivna samo okolica brez celice, v drugi samo celica brez okolice. Na robu med okolico in celico smo za robni pogoj določili gostoto tokaJmembrana po enačbi (3.5):
•               ty       _y       \ Imembrana                                                                                   (3.5)
J membrana        \    celica         okolica /7
"membrana
VceUca predstavlja električni potencial znotraj celice, VokoUca električni potencial v okolici celice, Ymembrana = 5×10-7 Sm-1 je specifična prevodnost membrane, d = 5 nm debelina membrane. V prvi fizikalni naravi, v kateri je aktiva le okolica, smo na robu nastavili +membrana, v drugi fizikalni naravi, v kateri je aktivna le celica, pa smo na istem robu nastavili - j membrana. S tem smo povezali VokoUca, ki je definiran v prvi fizikalni naravi, z VceUca, ki je definiran v drugi fizikalni naravi.
Za prevodnost znotrajcelične in izvencelične tekočine smo izbrali ytekocine = 0.2 Sm-1. Ker je robni pogoj funkcija spremenljivk, smo za reševanje sistema enačb uporabili nelinearni algoritem za izračun, in sicer smo uporabili iterativni nelinearni algoritem z SSOR (Symmetric Successive Over-Relaxation) metodo za začetno transformacijo (preconditioning) ter Quasi-Minimal Residual metodo za rešitev sistema enačb.
4 Rezultati
30
4 Rezultati
4.1   Uporaba skaliranih elektrod
V FEMLAB-u smo s pomočjo modela elektrod s premerom 2r = 0.7 mm ter dolžine 40 mm določili električno polje dveh igelnih elektrod z medsebojno oddaljenostjo d = 2 mm, 4 mm, 6 mm, 8 mm in 10 mm. Na elektrodi je bil priključen električni potencial V0 = ±100 V. Določili smo ekvipotencialne ploskve za Vs = ±25 V0, ±50 V0 ter ±75 V0. Rezultati so podani na sliki 4.1. Iz slike je razvidno, da sta ekvipotencialni ploskvi z Vs = ±25 V0, kar odgovarja 1 V0, izrazito deformirani zaradi nehomogenosti snovi. Zato bi bila napaka pri nadomestitvi takšne ekvipotencialne ploskve z valjem večja kot v primeru Vs = ±50 V, kar ustreza 1 V0. Seveda je napaka najmanjša v primeru Vs = ±75 V0, a hkrati je tudi faktor skaliranja k med nadomestno elektrodo ter originalno elektrodo najmanjši. Naš cilj je čim večji faktor skaliranja k ob sprejemljivi napaki nadomestitve ekvipotencialne ploskve z valjem, saj bo tako elektroda debelejša in gradnja mreže bolj enostavna. Iz primerjav slik lahko opazimo, da je, poleg faktorja skaliranja k, pomemben podatek o nadomestnih elektrodah tudi poznavanje razdalje med nadomestnima elektrodama ds ter razdalje med središčema nadomestnih elektrod d0s. Vidimo namreč lahko, da ekvipotencialna ploskev ne leži simetrično okrog elektrode, ampak je premaknjena stran od nasprotne elektrode. Za določanje nadomestnih elektrod je torej potrebno poznati tudi ta podatek. Podobno kot smo definirali faktor skaliranja k v enačbi (3.1), definiramo tudi faktor skaliranja razdalje med središčema elektrod l:
l = 100-( d—1).                                                                                              (4.1)
d0
Na podlagi slike 4.1 smo določili povprečni faktor skaliranja polmera k ter povprečni faktor skaliranja razdalje l. Rezultati so podani v tabeli (4.1). Zaradi prevelike razlike med eno in drugo ekvipotencialno ploskvijo za razmerje Vs = ±25 V0 ne moremo več določiti k in l. Na podlagi tega sklepamo, da je za izbrano razmerje prevodnosti n = 0.1 Sm-1 ter y2 = 1 Sm-1, sprejemljiva meja nadomestnih elektrod le do razmerja Vs / V0 = 0.5.
4 Rezultati
31
Slika 4.1: Ekvipotencialni ploskvi za Vs = ±75 V0 (3 električnega potenciala elektrod), ±50 V0 (1 električnega
potenciala elektrod) ter ±25 V0 (1 električnega potenciala elektrod) v okolici dveh igelnih elektrod sta prikazani
levo, v sredini ter desno. Elektrodi imata premer 2r = 0.7 mm, dolgi sta 40 mm, razdalja med njima d = 4 mm v
prvi vrstici, 6 mm v drugi, 8 mm v tretji ter 10 mm v četrti.
4 Rezultati
32
Iz slike 4.1 je razvidno, da se tudi pri razmerju Vs / V0 = 0.5 skaliranje elektrod za pozitivno ter za negativno elektrodo med seboj razlikujeta, saj se polmera leve in desne elektrode med seboj razlikujeta. Vzrok nesimetričnosti je v nehomogenosti tkiva. Če bi izbrali homogeno tkivo (?1 = ?2), bi bilo skaliranje obeh elektrod popolnoma enako.
Tabela 4.1: Faktor skaliranja polmera elektrod k ter faktor skaliranja razdalje l nadomestnih elektrod.
d [mm]	4		6		8		10	
Vs [V]	75	50	75	50	75	50	75	50
r [mm]	0.35	0.35	0.35	0.35	0.35	0.35	0.35	0.35
d0 [mm]	4.7	4.7	6.7	6.7	8.7	8.7	10.7	10.7
d0s [mm]	5.02	5.54	6.8	7.14	8.76	9.36	10.98	11.48
ds [mm]	3.81	2.96	5.37	3.85	7.1	5.75	9.14	6.99
rs [mm]	0.60	1.29	0.71	1.64	0.83	1.80	0.92	2.24
k = 100(rs / r-1)	72.86	268.57	104.28	370	137.14	415.71	162.86	541.43
l = 100(d0s / d0)	6.81	17.87	1.4925	6.57	0.69	7.59	2.62	7.29
Na podlagi rezultatov iz tabele 4.1 smo narisali grafa na sliki 4.2. Predstavljata nam faktor skaliranja polmera elektrod k ter razdalje med elektrodama l v odvisnosti od razmerja med Vs in V0. Za večjo preglednost smo dodali tudi izhodiščno točko, ko elektrod ne skaliramo in je vrednost obeh faktorjev v njej 0. Kakor je razvidno iz slike 4.2, so podatki na sliki podani kot razmerja, tako da veljajo za različne premere elektrod, razdalje in električne potenciale.
600 500 400	faktor skaliranja polmera k	
		
		
		¦r/d0=0.07
300		¦   r/d0=0.05
		A   r/d0=0.04
200		9   r/d0=0.03
		
0 1		
		
	0.75                                                  0.5        Vs/V	
faktor skaliranja raz dalje l
1
r/d0=0.07 r/d0=0.05 r/d0=0.04 r/d0=0.03
Slika 4.2: Faktor skaliranja polmera elektrod k (levo) ter faktor skaliranja razdalje med elektrodama l (desno) v
odvisnosti od izbranega razmerja med Vs ter V0.
Iz slike 4.2 je razvidno, da faktor skaliranja polmera elektrod k z manjšanjem razmerja med električnim potencialom ekvipotencialnih ploskev ter elektrod Vs / V0 izrazito narašča. Že za največje razmerje med polmerom in razdaljo med elektrodama r / ds, kjer je vpliv skaliranja
Vs/V
0.75
0.5
4 Rezultati
33
najmanjši, imamo v primeru Vs / V = 0.75 vrednost k ~ 100, pri vrednosti Vs / V0 = 0.5 pa že ~ 250. Poleg tega je izrazito tudi naraščanje k z manjšanjem razmerja med polmerom in razdaljo med elektrodama r / ds, drugače rečeno – k narašča z večanjem razdalje med elektrodama. Obratno je zaznati za faktor l, saj njegova vrednost z večanjem razdalje med elektrodama izrazito upada in pri razmerju r / ds = 0.05 ne preseže vrednosti 8.
Slika 4.3: Potek električnega potenciala V v y osi izhodiščnega modela je prikazan s polno črto, s prekinjeno črto
pa je prikazan potek električnega potenciala v osi y modela s skaliranimi elektrodami. Os y je os, ki poteka skozi
obe elektrodi na polovici dolžine elektrod. Razdalja med elektrodama je 10 mm, razmerje Vs / V0 = 0.75. Na levi
strani je prikazan potek potenciala v modelu na slikah 4.1 do 4.4, v katerem sta snovi z dvema različnima
prevodnostma, na desni pa je predstavljen primer za homogen model.
Za Vs / V = 0.75 smo izračunali model s skaliranimi elektrodami in rezultate primerjali z rezultati izhodiščnega modela. Primerjava rezultatov je predstavljena na sliki 4.3. Prikazan je potek električnega potenciala na osi y. Os y je tista os, ki poteka skozi obe elektrodi na polovici dolžine elektrod. S polno črto je prikazan potek električnega potenciala v izhodiščnem modelu, s prekinjeno črto pa je prikazan potek električnega potenciala v modelu s skaliranimi elektrodami. Iz slike 4.3 levo, kjer je predstavljen potek potenciala v modelu z enakimi lastnostmi snovi kot na sliki 4.1, lahko vidimo, da je napaka pri nehomogenem tkivu opazna. V sami bližini elektrod je napaka manjša, bolj izrazita pa na robu. Vzrok je v tem, da smo tako za levo kot tudi za desno elektrodo uporabili isti faktor skaliranja, kar je le približek. Velikost volumna, ki ga objame ekvipotencialna ploskev v predelu modela, kjer je prevodnost višja, je večja kot v drugem delu. V homogenem modelu na desni sliki pa se rezultati le malo razlikujejo.
Ker je za ocenjevanje učinkov na organizme bolj kot električni potencial, predstavljen na sliki 4.3, pomembna električna poljska jakost, je absolutna vrednost le-te prikazana na sliki 4.4.
4 Rezultati
34
Slika 4.4: Enako kot na sliki 4.3, le da je pokazan potek absolutne vrednosti električne poljske jakosti E
Iz slike 4.4 je razvidno, da je napaka električne poljske jakosti tudi v večji oddaljenosti od izhodišča precej manjša kot napaka električnega potenciala. To se še bolj natančno vidi na sliki 4.5, ki prikazuje povečano sliko 4.4 levo.
Slika 4.5: Povečana slika 4.4 levo. Očitno je, da je napaka v porazdelitvi absolutne vrednosti električne poljske
jakosti E med izhodiščnim modelom (polna črta) ter modelom s skaliranimi elektrodami (prekinjena črta),
majhna.
4 Rezultati
35
4.2 Model možganov s tumorjem
Iz CT slik smo zgradili model človeških možganov s tumorjem. V model smo vključili tudi vzporedne osemigelne elektrode, kakršne se uporabljajo pri Cliniporatorju. Model smo brez uporabe nelokalnega sklapljanja rešili le, če smo elektrode iz modela izrezali. Zaradi majhnega premera elektrod je namreč gradnja mreže v elektrodah zelo zahtevna, saj so okrogle in imajo majhen premer v primerjavi z velikostjo celega modela. Poleg tega nas porazdelitev električnega polja znotraj elektrod ne zanima, ker je glede na njihovo veliko električno prevodnost električna poljska jakost v njih približno nič. Porazdelitev električne poljske jakosti v prerezni ravnini x = 55 mm je predstavljena na sliki 4.6.
Slika 4.6: Porazdelitev električne poljske jakosti v prerezni ravnini x = 55 mm za celoten model možganov s tumorjem. Ker so elektrode izrezane, znotraj njih ni polja in so zato bele. Barvna skala na desni strani je v V/m.
Model smo rešili tudi z nelokalnim sklapljanjem. Pri tem elektrod ni bilo potrebno izrezati iz modela, saj gradnja mreže tako v večji kot tudi v manjši geometriji ni predstavljala večjih težav. Porazdelitev električne poljske jakosti v prerezni ravnini x = 55 mm je predstavljena na sliki 4.7. Za boljšo primerjavo je na sliki 4.8 predstavljena porazdelitev električne poljske jakosti in električnega potenciala samo za področje manjše geometrije.
4 Rezultati
36
Slika 4.7: Porazdelitev električne poljske jakosti v prerezni ravnini x = 55 mm za model možganov s tumorjem in nelokalnim sklapljanjem. Zgornja slika predstavlja veliko geometrijo, spodnja pa malo, izrezano geometrijo.
Barvna skala na desni strani je v V/m.
4 Rezultati
37
Slika 4.8: Porazdelitev električne poljske jakosti (zgoraj) in električnega potenciala (spodaj) v prerezni ravnini x
= 55 mm za model možganov s tumorjem in nelokalnim sklapljanjem. Levi sliki predstavljata porazdelitev električne poljske jakosti izhodiščnega modela, desni dve pa v modelu z nelokalnim sklapljanjem. Barvna skala
na desni strani je v V/m.
Iz slike 4.8 je razvidno, da se tako porazdelitev električne poljske jakosti kakor tudi električnega potenciala v bližini elektrod ujema. Zaradi bolj natančne mreže je rezultat za porazdelitev električne poljske jakosti na desni sliki, ki predstavlja model z nelokalnim sklapljanjem, bolj podroben kot na levi.
4 Rezultati
38
Slika 4.9: Porazdelitev električnega potenciala v prerezni ravnini x = 140 mm za izhodiščni model (zgoraj) in model z nelokalnim sklapljanjem (spodaj). Barvna skala na desni strani je v V.
Poleg porazdelitve električne poljske jakosti in električnega potenciala v bližini elektrod je na sliki 4.9 predstavljena porazdelitev električnega potenciala v prerezni ravnini x = 140 mm, kjer, kakor v bližini elektrod, ni opaziti večjih razlik med modeloma.
4 Rezultati
39
4.3 Model celice
Uporabo nelokalnega sklapljanja smo preizkusili tudi na modelu celice. Osnova geometrije je krogla polmera 10 µm v središču kocke z robom 40 µm. Krogla predstavlja celico, kocka pa okolico. Na obeh skrajnih koncih celice glede na os y, kjer so pričakovane največje vrednosti električne poljske jakosti, smo izrezali dva manjša dela celice osnovne geometrije, velika 100 nm × 30 nm × 100 nm. Oba izrezana dela smo nadomestili z geometrijama plus in minus, sestavljenima iz treh kvadrov: 100 nm × 30 nm × 100 nm (okolica), 100 nm × 20 nm × 100 nm (celična membrana) ter 100 nm × 15 nm × 100 nm (notranjost celice). Na meji rezanja smo s pomočjo spremenljivke sklapljanja določili sklopljen robni pogoj. Porazdelitev električnega potenciala za primer tokovnega polja je predstavljena na slikah 4.10 in 4.11.
Slika 4.10: Porazdelitev električnega potenciala v modelu celice. Prikazan je potek električnega potenciala v
ravnini x = 0. Na desnem robu je priključen električni potencial 1 V, na levem električni potencial -1 V. Okolica
drobnih pik na skrajnih koncih celice glede na os y, ki sta izrezana dela, je predstavljena na sliki 4.11 zgoraj,
sama izrezana dela na sliki 4.11 spodaj. Barvna skala na desni strani je v V.
4 Rezultati
40
Slika 4.11: Porazdelitev električnega potenciala v modelu celice v ravnini x = 0. Zgornji dve sliki predstavljata povečana dela slike 4.10 v okolici izrezanih delov. Na spodnjih dveh slikah je predstavljen električni potencial v
izrezanih delih. Barvna skala na desni strani je v V.
4 Rezultati
41
Iz slike 4.11 lahko vidimo, da je električni potencial zaradi sklopljenih robnih pogojev zvezen. Še bolje je to razvidno iz slike 4.12, kjer smo spodnji dve porazdelitvi na sliki 4.11 vstavili v zgornji dve in dobili sestavljeno porazdelitev električnega potenciala iz velike geometrije, geometrije minus in geometrije plus.
Slika 4.12: Porazdelitev električnega potenciala v modelu celice v ravnini x = 0. Sliki sta sestavljeni iz spodnjih dveh ter zgornjih dveh slik slike 4.11. Barvna skala na desni strani je v V.
Napaka v porazdelitvi potenciala je opazna le na zgornjem in spodnjem robu, kjer se stikata izrezani del (geometrija minus ali plus) ter velika geometrija. Ker z razliko od izrezanih delov v veliki geometriji ni membrane, je bilo potrebno na prerezu membrane namesto sklapljanja definirati izolativne pogoje, kar je povzročilo napako na tem robu. Zanimivo je, da kljub sorazmerno veliki prevodnosti znotrajcelične tekočine, ?tekocine = 0.2 Sm-1 v primerjavi s prevodnostjo membrane ?membrana = 5×10-7 Sm-1, lahko iz slike zanamo razliko v električnem potencialu znotraj celice v geometiji plus in geometriji minus, a je ta razlika majhna ter posledično tudi električna poljska jakost znotraj celice zelo majhna.
5 Razprava
42
5 Razprava
V okviru naloge smo preučili možnost uporabe dveh načinov, s katerima bi lahko rešili problem velikih geometrijskih razmerij pri modeliranju in izračunu z numeričnimi metodami in jih do sedaj v Laboratoriju za biokibernetiko še nismo uporabljali.
Prvi način je skaliranje elektrod. Pri vseh problemih, v katerih kot izvor električnega polja služijo elektrode, nameščene v tkivu oz. na tkivu, je gradnja mreže ponavadi zahtevna. Elektrode so namreč zelo tanke (premer manj kot 1 mm) v primerjavi z velikostjo modela. V nalogi smo z modeli ocenili, za kolikšen faktor skaliranja lahko v modelu povečamo elektrode ter temu ustrezno znižamo napetost na njih, da bo rezultat še dovolj dober.
Drugi način je uporaba nelokalnega sklapljanja (non-local coupling). Pri metodi končnih elementov ali končnih diferenc na podlagi zgrajene mreže sestavimo matrično enačbo (1.1), ki jo je potrebno rešiti. Namesto da bi imeli samo eno geometrijo, smo iz celotne geometrije izrezali drugo, manjšo, ki zavzema tisti del modela, ki predstavlja največ težav za gradnjo modela oziroma smo v njem želeli rešitev z večjo prostorsko ločljivostjo. Robne pogoje na mestu rezanja smo določili tako, da so se preslikali iz večje geometrije na manjšo in obratno. Tako smo en model sestavili iz dveh geometrij, matrična enačba pa je bila zaradi preslikovanja robnih pogojev ena sama.
5.1 Uporaba skaliranih elektrod
Metodo skaliranja elektrod smo uporabili na paru igelnih elektrod premera 0.6 mm in dolžine 40 mm. Rezultati v obliki ekvipotencialnih ploskev z električnim potencialom Vs = 0.25, 0.5 in 0.75 V0 so predstavljeni na sliki 4.1. V nehomogenem modelu z razmerjem med prevodnostma ?2 / ?1 = 10 je opaziti nesimetričnost ekvipotencialnih ploskev. Namreč, ekvipotencialne ploskve okrog elektrode, ki je v večjem delu v snovi z višjo prevodnostjo, objamejo večji volumen kot tiste okrog elektrode, ki je v večjem delu v snovi z nižjo prevodnostjo. Pri določanju nadomestnih skaliranih elektrod smo uporabili povprečni premer obeh ekvipotencialnih ploskev. Zaradi tega se odziv na sliki 4.3 levo, kjer je prikazan potek
5 Razprava
43
električnega potenciala V po osi y, razlikuje. Pomembnejše je dejstvo, da se potek absolutne vrednosti električne poljske jakosti E, ki je predstavljen na sliki 4.4, zelo malo razlikuje. Učinke električnega polja na žive organizme, tkiva in celice namreč bolj kot sam električni potencial določa absolutna vrednost električne poljske jakosti, saj absolutna vrednost električne poljske jakosti na različnih strukturah, kot denimo celični membrani, povzroči vsiljeno napetost in s tem učinke na celice ali biološka tkiva. Porazdelitev absolutne vrednosti električne poljske jakosti se v modelu s skaliranimi elektrodami ujema s porazdelitvijo v izhodiščnem modelu.
V  splošnem bi bilo potrebno za natančen izračun porazdelitve električne poljske jakosti in električnega potenciala v modelu izračun opraviti v neskončnem prostoru. Ker to seveda ni mogoče, model, v katerem imamo lokalni vir, zapremo v dovolj veliko območje in na robu tega območja lahko z majhno napako predpostavimo izolativni robni pogoj. Večje kot je območje, manjša bo napaka zaradi omejitve območja računanja. V tistih primerih, ko nas zanima samo bližnje polje ob elektrodah, skaliranje ni smiselno. Bližnje polje je v tem primeru mišljeno kot električno polje, eno ali dve razdalji med elektrodama oddaljeno od elektrod. Čeprav moramo v primeru, ko ne uporabimo skaliranja elektrod, uporabiti manjše območje, to na porazdelitev električnega polja v bližini elektrod ne vpliva veliko. Če nas zanima porazdelitev električnega polja tudi v večji oddaljenosti od elektrod, je skaliranje smiselno. S skaliranjem elektrod lahko sorazmerno s povečanjem polmera elektrod povečamo tudi velikost območja in s tem dosežemo manjšo napako zaradi omejitve območja računanja.
V  modelu smo skaliranje elektrod preizkusili za razmerje v prevodnosti y2 / n = 10. V splošnem v tkivih običajno ne srečamo večjih razmerij v prevodnosti, razen v primeru kože. Koža ima nekaj velikostnih razredov manjšo prevodnost kot ostala tkiva. Ker pa je koža zelo tanka, jo igle v primeru, da so zapičene v tkivo, prebadajo in zato na porazdelitev električne poljske jakosti in električnega potenciala v samem tkivu koža ne vpliva. Na zunanji plasti kože je namreč vedno izolativni robni pogoj, ki se zaradi zelo majhne prevodnosti kože pravzaprav prenese tudi na notranjo plast kože. Poleg tega elektrode prebadajo kožo pravokotno na elektrode. Če so namreč elektrode pravokotne na mejo dveh prevodnikov in jo prebadajo, je vpliv različnih prevodnosti zelo majhen, še več, če so elektrode ravno polovico v enem prevodniku in ravno polovico v drugem prevodniku, se zaradi simetričnosti vpliv različnih prevodnosti za izotropno tkivo izniči. Previdnost pri uporabi skaliranja elektrod je potrebna v primeru, da poteka meja med dvema prevodnikoma med elektrodama ali je večina prve elektroda v enem prevodniku in večina druge elektrode v drugem prevodniku. V tem
5 Razprava
44
primeru je razlika v volumnu, ki ga objamejo ekvipotencialne ploskve okrog prve elektrode ter ekvipotencialne ploskve okrog druge elektrode največja ter posledično tudi največja napaka zaradi uporabe skaliranih elektrod.
Iz podatkov na sliki 4.2 lahko za par elektrod z danim polmerom ter medsebojno razdaljo razberemo faktor skaliranja polmera elektrod k ter faktor skaliranja razdalje med elektrodama l za razmerja med prevodnostma ?2 / ?1 ? 10. Uporabne vrednosti faktorja k , torej tiste vrednosti, kjer je napaka zaradi skaliranja dovolj majhna, se gibljejo med 300 in 500. Tako lahko z metodo skaliranja elektod povečamo premer elektrod v modelu 3 do 5 krat. Za večja razmerja v prevodnosti ter za drugačno geometrijo elektrod bi bilo potrebno pred uporabo metode z modelom ugotoviti, kakšna geometrija elektrod bi ustrezno nadomestila ekvipotencialne ploskve okrog elektrod.
5.2 Model možganov s tumorjem in model celice
Metodo nelokalnega sklapljanja smo preizkusili na modelu človeških možganov s tumorjem in vzporednimi 8 igelnimi elektrodami, kakršne se uporabljajo pri Cliniporatorju. Rešitev modela z nelokalnim sklapljanjem smo primerjali z rešitvijo izhodiščnega modela. Iz slik 4.6 do 4.8 je razvidno, da se tako rezultati za porazdelitev električne poljske jakosti kakor tudi električnega potenciala med obema modeloma dobro ujemajo.
Drugi model, kjer smo uporabili nelokalno sklapljanje, je model celice z membrano. Sicer bi bilo mogoče zgraditi model celice s skalirano membrano, tako da bi njeno debelino povečali za nek dovolj velik faktor (denimo 200) in temu ustrezno priredili prevodnost membrane, vendar smo želeli preveriti, ali metoda nelokalnega sklapljanja omogoča izgradnjo modela celice z realno debelino membrane. Pri gradnji modela smo se srečali s številnimi težavami. Ker je v FEMLAB-u mogoče nelokalno sklapljanje določiti le na mejah rezanja, ki so vzporedne z ravninami xy, xz in yz, smo lahko v drugo geometrijo izrezali le en manjši del celične membrane, in sicer smo se odločili tak izrez narediti na dveh kapah celice, kjer je električna poljska jakost največja. Izbrati smo morali zelo majhen izrez, saj smo v izrezani geometriji morali izrezane dele nadomestiti s kvadri, sicer gradnja mreže ni bila mogoča (glej sliki 3.5 in 3.6). V primeru velikega izrezanega dela bi bila napaka zaradi uporabe nadomestne geometrije v obliki kvadrov namesto krogelnih odsekov prevelika. Tudi z izračunom so bile številne težave, saj je bilo zaradi odvisnosti robnih pogojev od vrednosti
5 Razprava
45
iskanih spremenljivk potrebno za reševanje sistema enačb uporabiti nelinearni algoritem. Zahtevnost reševanja sistema enačb določa število prostostnih stopenj modela (degrees of freedom ali DOF), ki pove, koliko neznanih spremenljivk vsebuje sistem enačb. Število prostostnih stopenj modela je odvisno predvsem od števila vozlišč v modelu in stopnje elementov, uporabljenih v modelu. Stopnja elementov je stopnja polinoma, s katerim opisujemo potek posamezne spremenljivke v posameznem elementu. Kljub majhnemu številu prostostnih stopenj modela, 27740, je bilo potrebno za rešitev uporabiti iterativni algoritem, saj direktni algoritem ni deloval, a tudi iterativni algoritem je imel precej težav s konvergenco.
Glede na dejstvo, da se rezultati za model možganov s tumorjem brez nelokalnega sklapljanja in z nelokalnim sklapljanjem (slike 4.6 do 4.9) skorajda ne razlikujejo ter glede na številne težave pri modelu celice z membrano, se pojavlja vprašanje, v čem je smiselnost uporabe nelokalnega sklapljanja, ki zahteva dodatno pozornost ter delo. Odgovor se skriva v prednostih metode nelokalnega sklapljanja:
•    v modelu brez uporabe nelokalnega sklapljanja smo zaradi majhnega premera elektrod in posledično zahtevne gradnje mreže v elektrodah lahko mrežo zgradili le, če smo elektrode iz modela izrezali, v modelu z nelokalnim sklapljanjem pa smo mrežo lahko zgradili tudi v elektrodah;
•    model brez nelokalnega sklapljanja je imel število prostostnih stopenj 61262. Ker gostejše mreže, kot smo jo zgradili, imela je 8916 vozlišč, zaradi zahtevnosti geometrije ni bilo mogoče zgradi, smo v njem uporabili kvadratne elemente, saj ti ob istem številu vozlišč dajo bolj natančen rezultat kot linearni elementi. Model z nelokalnim sklapljanjem je imel 5559 vozlišč v veliki geometriji in 72286 vozlišč v mali geometriji, skupno skoraj 10 krat več kot model brez nelokalnega sklapljanja. Število prostostnih stopenj je bilo 77845. Zaradi bolj podrobne mreže v primeru nelokalnega sklapljanja v okolici elektrod je rezultat za porazdelitev električne poljske jakosti bolj podroben kot v modelu brez uporabe nelokalnega sklapljanja, kar je razvidno tudi iz slike 4.8 zgoraj;
•    zaradi zelo zahtevne mreže je v modelu brez nelokalnega sklapljanja zelo težko premikati elektrode. Da se je mreža sploh lahko zgradila, je bilo potrebno na sredini med elektrode dodati prerezno ravnino in nato zgraditi mrežo. Položaj prerezne ravnine je bilo potrebno določiti s preizkušanjem, kdaj bo avtomatskemu algoritmu za gradnjo mreže le-to uspelo zgraditi. V modelu z nelokalnim sklapljanjem je bila
5 Razprava
46
izgradnja mreže dosti manj problematična, premikanje elektrod pa je bilo še vedno zelo omejeno. Nastajal je namreč problem z analizo geometrije, kar je posledica dejstva, da je geometrija zgrajena iz slik. V primeru, da bi bila geometrija sestavljena iz enostavnih geometrijskih teles, tega problema ne bi bilo;
•    kljub številnim težavam z geometrijo v modelu celice z membrano je bila velikost izrezanega dela še vedno dovolj velika, da bi v njem lahko modelirali poro. Problem računanja modela s poro bi gotovo bil v konvergenci izračuna, ki je bila že sedaj zelo počasna.
Slabosti nelokalnega sklapljanja so predvsem v tem, da:
•    morajo biti v FEMLABU meje rezanja, na katerih določimo sklapljanje, vzporedne ravninam xy, xz in yz, lahko pa je koordinatni sistem poljubno postaljven;
•    se prostostne stopnje v geometrijah seštevajo. V primeru modela možganov s tumorjem z nelokalnim sklapljanjem to pomeni, da je bilo število prostostnih stopenj modela 5559 +72286 = 77845;
•    je pri običajnih problemih brez nelokalnega sklapljanja matrika za reševanje redka. Vrednost iskane spremenljivke v posameznem vozlišču je namreč ponavadi odvisna le od vrednosti v nekaj sosednjih vozliščih, vsi ostali elementi v vrstici pa so enaki nič. Takšna matrika je za reševanje precej bolj enostavna ter zahteva manj računskega časa kot manj redka matrika, ki nastane pri uporabi nelokalnega sklapljanja. Včasih je razlika zelo majhna, lahko pa je tudi zelo velika. FEMLAB ima za ta namen sicer vgrajen poseben algoritem za neredke matrike, ki deloma odpravi težavo;
•    je pri uporabi nelokalnega sklapljanja dosti večja verjetnost, da bo imel problem slabšo konvergenco;
•    je trenutno edini nam znani programski paket, ki ima predvideno nelokalno sklapljanje, FEMLAB.
6 Zaključek
47
6 Zaključek
Razvoj programskih paketov za modeliranje in izračun z numeričnimi metodami ter hitro večanje računalniških zmogljivosti omogoča vedno bolj natančno računanje. Kar je danes ravno še mogoče, bo že v bližnji prihodnosti nekaj vsakdanjega. A s tem metode, ki rešujejo problem velikih geometrijskih razmerij, ne bodo pozabljene, saj se bomo pri raziskovalnem delu vedno srečevali s problemi, kjer bodo prisotna velika geometrijska razmerja ali bomo želeli večjo prostorsko ločljivost rešitve.
Metoda skaliranja elektrod je preprosta metoda, a žal z omejenimi možnostmi. Uporabna je le v modelih, v katerih uporabljamo tanke elektrode v primerjavi z velikostjo modela. Z njo lahko razmerje med premerom elektrod in velikostjo modela zmanjšamo štiri do pet krat. Sam faktor skaliranja ni velik, a vendar pomeni, da lahko namesto modela, velikega npr 100 ×100 × 100 mm z elektrodami izračunamo model velik 400 ×400 × 400 mm z nadomestnimi elektrodami. V veliki večini primerov, ko računamo porazdelitve električnega polja v okolici elektrod, ki se uporabljajo pri elektroporaciji in genski transfekciji, je takšna velikost modela dovolj velika. Metodo omejuje dejstvo, da je potrebno za različno število in postavitev elektrod določiti nove nadomestne elektrode.
Metoda nelokalnega sklapljanja (non-local coupling) je v osnovi preprosta metoda, ki pa jo mora podpirati programski paket, v katerem model zgradimo in računamo. Žal je to po nam znanih podatkih trenutno le FEMLAB. Je zelo splošna, zato je uporabna na različnih področjih, tako v biomedicinski tehniki kakor tudi v tehniki na sploh. Uporabna je tako v geometrijah z elektrodami, kakor metoda skaliranja elektrod, kakor tudi v vseh drugih geometrijah, kjer se srečujemo z velikim razmerjem med velikostjo celotne geometrije in velikostjo njenih sestavnih delov. Poleg tega metoda nelokalnega sklapljanja omogoča v nekem področju geometrije izračun z večjo prostorsko ločljivostjo kot v celotnem modelu ter določanje periodičnih robnih pogojev. V FEMLAB-u je omejena je z dejstvom, da morajo biti ravnine rezanja vzporedne s koordinatnim sistemom, ki pa je lahko poljubno postavljen, po drugi strani pa nobeden od drugih programskih paketov za modeliranje in izračun z numeričnimi metodami nelokalnega sklapljanja niti ne omogoča.
7 Literatura
48
7 Literatura
Biotehniška fakulteta, Univerza v Ljubljani. Študij biotehnologije. Biotehniška fakulteta, Univerza v Ljubljani, 2004. http://www.bf.uni-lj.si/stud_prg/biotehnologija.pdf. (1. 6. 2004)
Comsol AB. FEMLAB Reference Manual v. 2.2. Comsol AB., Stockholm, 2001.
Kotnik T, Miklavčič D, Analytical description of transmembrane voltage induced by electric fields on spheroidal cells. Biophys. J. 79: 670-679, 2000.
Pavlin M, Miklavčič D, Biophysical Effective Conductivity of a Suspension of Permeabilized Cells: A Theoretical Analysis. Biophys. J. 85: 719-729, 2003.
Pavšelj N, Porazdelitev električne poljske jakosti v podkoćnem tumorju pred in po elektroproaciji celic. Magistrska naloga, Fakulteta za elektrotehniko, Univerza v Ljubljani, 2002.
Sinigoj AR, ELMG polje. Založba FE, Ljubljana, 1996.
Susil R, Šemrov D, Miklavčič D, Electric field induced transmembrane potential depends on cell density and organization. Electro. Magnetobiol. 17: 391-399, 1998.
Šel D, Določitev časovno diskretnega modela elektropermeabilizacije tkiva. Doktorska disertacija, Fakulteta za elektrotehniko, Univerza v Ljubljani, 2003.
Sel D, Mazeres S, Teissié J, Miklavcic D. Finite-element modeling of needle electrodes in tissue from the perspective of frequent model computation. IEEE Trans. Biomed. Eng. 50: 1221-1232, 2003.
Univerza v Ljubljani. Biomedicina Univerzitetni znanstveni podiplomski študij. Univerza v Ljubljani, Ljubljana, 2003. http://www.uni-lj.si/Biomedicina/Uvod.pdf. (1. 6. 2004)
7 Literatura
49
Valič B, Golzio M, Pavlin M, Schatz A, Faurie C, Gabriel B, Teissié J, Rols MP, Miklavčič D. Effect of electric field induced transmembrane potential on spheroidal cells: theory and experiment. Eur. Biophys. J. 32: 519-528, 2003.
The Whitaker Foundation. Definition of Biomedical Engineering. The Whitaker Foundation, 2004. http://www.whitaker.org/glance/definition.html (1. 6. 2004)
Zupančič B, Karba R, Matko D. Simulacije dinamičnih sistemov. Založba FE, Ljubljana, 1995.
50
Izjava
Izjavljam, da sem magistrsko nalogo izdelal samostojno pod vodstvom mentorja prof. dr. Damijana Miklavčiča, dipl. ing. Izkazano pomoč drugih sodelavcev sem v celoti navedel v zahvali.
Blaž Valič