Ontwikkeling van numerieke methoden voor thermotherapie van leverkanker

Ontwikkeling van numerieke methoden voor thermotherapie van leverkanker Frederik Soetaert

Promotoren: prof. dr. ir. Luc Dupré, dr. ir. Guillaume Crevecoeur Masterproef ingediend tot het behalen van de academische graad van Master in de ingenieurswetenschappen: toegepaste natuurkunde

Vakgroep Elektrische energie, Systemen en Automatisering Voorzitter: prof. dr. ir. Jan Melkebeek Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2011-2012

Ontwikkeling van numerieke methoden voor thermotherapie van leverkanker Frederik Soetaert

Promotoren: prof. dr. ir. Luc Dupré, dr. ir. Guillaume Crevecoeur Masterproef ingediend tot het behalen van de academische graad van Master in de ingenieurswetenschappen: toegepaste natuurkunde

Vakgroep Elektrische energie, Systemen en Automatisering Voorzitter: prof. dr. ir. Jan Melkebeek Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2011-2012

De auteur geeft de toelating deze masterproef voor consultatie beschikbaar te stellen en delen van de masterproef te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze masterproef. The author gives permission to make this master dissertation available for consultation and to copy parts of this master dissertation for personal use. In the case of any other use, the limitations of the copyright have to be respected, in particular with regard to the obligation to state expressly the source when quoting results from this master dissertation. Gent, 1 juni 2012 De auteur

Frederik Soetaert

Voorwoord “Burgerlijk ingenieur in de kernfysica,” luidde mijn antwoord toen ik zes jaar oud was op de vraag wat ik later wou worden. Wist ik veel dat deze studierichting niet eens bestond. Mijn antwoord was enerzijds gebaseerd op het aureool van respect dat op volwassenen hun gezicht verscheen als ze het woord “burgerlijk ingenieur” in de mond namen. Het kernfysica-aspect werd anderzijds ingegeven door de lezing van een krantenartikel over het fabriceren van een nieuw chemisch element en de zoektocht naar het stabiele eiland der superactiniden. Ondanks mijn uiterst elementaire kennis op dat moment, was ik blijkbaar toen al enorm ge¨ınteresseerd in de minst populaire deelgebieden van menselijke activiteiten. Het was dan ook niet verwonderlijk dat ik een opleiding tot burgerlijk ingenieur in de toegepaste natuurkunde aanvatte. De meest fascinerende, complexe en fundamentele aspecten van de natuur zijn de laatste vijf jaar op mijn radar verschenen. Dat deze periode ook gepaard ging met een intensief studieschema nemen we er retrospectief graag bij. Om prof. Waroquier te citeren bij aanvang van het examen Kwantummechanica II: “Onthou één ding: het is de mooiste, maar ook de moeilijkste opleiding.” De essentie van vijf jaar in één zin gebald. Ik ben dankbaar dat de UGent als enige in Vlaanderen deze opleiding aanbiedt. Het verandert je denkbeeld ingrijpend en laat je terugvallen op een stevige fundering die je in de meest uiteenlopende studiegebieden kan aanwenden. Net wat een universitaire opleiding zou moeten beogen: een brede basiskennis die je toelaat om complexe materie uit allerlei invalshoeken zelfstandig te kunnen bestuderen. Sapere aude. Naast mijn voorliefde voor ingenieurswetenschappen had ik steeds een grote interesse voor geneeskundige onderwerpen. Ik was dan ook heel tevreden dat ik het onderwerp van deze masterproef terugvond in de extensieve lijst van mogelijke masterbeproevingen. Dit onderwerp is werkelijk een overlap van verschillende van mijn interessegebieden. In de eerste plaats wil ik dan ook prof. dr. ir. Luc Dupré bedanken om mij de kans aan te bieden om me te verdiepen in dit mooie onderwerp. Ik ben hem eveneens dankbaar voor het vlotte verloop van mijn thesis. Ik ben verheugd deel uit te maken van een dergelijke kwalitatief hoogstaande en aangename vakgroep. Verder wens ik ook prof. dr. ir. Patrick Segers te bedanken om deel uit te maken van mijn beoordelingscommissie. Een speciale laudatio gaat uit naar mijn copromotor en dagelijkse begeleider, dr. ir. Guillaume Crevecoeur. Ik heb zelden het genoegen gehad een persoon te mogen ontmoeten die een dergelijk niveau van enthousiasme en kennis op dagelijkse basis aan boord weet te leggen in alles wat hij doet. Dankzij zijn input is dit werk uitgegroeid tot iets waar ik met enige trots kan op terugkijken.

iv

v Ik ben hem enorm dankbaar dat hij mijn kennis op tal van deelgebieden heeft vergroot: het programmeren in C++, het kritisch bestuderen van numerieke methoden, hoe objectieve analyses van numerieke simulaties te maken en hoe deze resultaten op een gestructureerde en duidelijke manier weer te geven. Ik kijk dan ook vol vertrouwen uit naar de komende jaren waarin ik mij verder mag ontwikkelen in deze materie onder zijn toeziend oog. De universiteit heeft me naast een opleiding ook een bepaalde omgeving aangereikt. Hier heb ik enkele bijzondere mensen mogen ontmoeten. Ik apprecieer de vele wetenschappelijke discussies, maar bovenal denk ik met plezier terug aan de vele schijnbare faits divers. Ondanks het feit dat deze periode in ons leven wordt afgesloten, ben ik zeker dat deze vriendschapsbanden tot ver na onze studies zullen doorlopen. Verder moet ik ook mijn dank betuigen aan andere vrienden en familieleden die me mooie momenten hebben bezorgd en in de toekomst nog zullen bezorgen. Ten slotte wens ik mijn ouders, broer en zus te bedanken. Per definitie zijn mijn ouders verantwoordelijk voor mijn bestaan. Hun invloed op mijn ontwikkeling strekt zich echter veel verder uit. Ook al besef ik het misschien niet altijd, ze doen ongelofelijk veel inspanningen om mijn leven gemakkelijker te maken. Zo hebben ze me van kindsbeen af steeds gesteund in mijn persistente zoektocht naar kennis en hebben ze zich moeite noch troost gespaard in de ondersteuning van mijn studies. Het was een plezier om op te groeien in een dergelijke omgeving: een ideale voedingsbodem voor geluk. Veel meer kan en mag een kind niet vragen van zijn ouders. Waarvoor mijn oprechte dank! Frederik Soetaert Gent, 1 juni 2012

Overzicht Ontwikkeling van numerieke methoden voor thermotherapie van leverkanker Frederik Soetaert Promotoren: prof. dr. ir. Luc Dupré, dr. ir. Guillaume Crevecoeur Begeleider: dr. ir. Guillaume Crevecoeur Masterproef ingediend tot het behalen van de academische graad van Master in de ingenieurswetenschappen: toegepaste natuurkunde Vakgroep Elektrische energie, Systemen en Automatisering Voorzitter: prof. dr. ir. Jan Melkebeek Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2011-2012

Samenvatting Traditionele therapieën zoals chemo- en radiotherapie zijn niet altijd optimaal voor het behandelen van leverkanker. Getuige daarvan zijn de hoge mortaliteitscijfers bij patiënten met leverkanker. Chirurgie is slechts bij een beperkt aantal patiënten een effectieve behandelstrategie. De medische wereld gebruikt thans in toenemende mate alternatieve strategieën zoals thermotherapie waarbij de tumoren met warmte bestreden worden. Bipolaire radiofrequente ablatie is een mogelijke vorm van thermotherapie en kan het leverweefsel lokaal opwarmen door elektrische stromen te injecteren met behulp van twee elektroden die door een stroomkring aangestuurd worden. Om deze behandelingsmethode beter te begrijpen dienen numerieke methoden ontwikkeld te worden. Het wiskundig model bestaat uit drie gekoppelde deelproblemen nl. het elektromagnetisch deelprobleem, het thermisch deelprobleem en de thermische kostenfunctie. Dit proefschrift ontwikkelt numerieke oplossingsmethoden voor radiofrequente ablatie waarbij een stabiliteitsanalyse uitgevoerd wordt, de invloed van (materiaal)parameters in het wiskundig model op de spatio-temporele temperatuursdistributie wordt nagegaan en hoe de huidige thermotherapie verbeterd kan worden. Numerieke simulaties tonen immers aan dat de standaard bipolaire radiofrequente ablatie een suboptimale behandelingsmethode is. We stellen een meer effectieve methode voor waarbij een gepulste stroomvorm gebruikt wordt waardoor een opbouw van temperatuur mogelijk is in de zone tussen de twee naalden. De ontwikkelde numerieke methoden geven aan dat het gebruik van een dergelijke therapie kan zorgen voor een verbeterde lokale ablatie van leverkanker. Trefwoorden: kankertherapie, numerieke methoden, stabiliteit, gekoppeld probleem, radiofrequente ablatie vi

Development of numerical methods for hyperthermia of liver cancer Frederik Soetaert Supervisor(s): Luc Dupré, Guillaume Crevecoeur Abstract—Liver cancer is one of the most difficult cancers to treat resulting in a very high mortality rate. Traditional treatments such as surgery, chemo- and radiotherapy are not always able to cure these cancers. An alternative approach, so-called radiofrequency ablation (RFA), applies an alternating current in the tumor that leads to a local heating of the tumoral tissue and thus removal of the cancer. Accurate numerical methods are needed so to understand these phenomena to a better extent. The work performed during the thesis presents numerical methods for simulating RFA and allows to study the impact of material properties on the spatio-temporal temperature distribution. Furthermore, we propose a more optimal RFA treatment using pulsed currents that could potentially increase the efficacy of RFA in clinical practice. Keywords—cancer therapy, numerical methods, stability, coupled problems, radiofrequency ablation

I. I NTRODUCTION RADITIONAL therapies for cancer are mostly a combination of surgery, chemo- and radiotherapy. These cancer therapies are not always an adequate or effective treatment. This occurs especially when treating liver cancers, one of the most occurring forms of cancer. Primary liver tumors, so-called hepatocellular carcinoma (HCC), occur approximately in 740000 cases per year. Unfortunately, the five year survival rate for HCC is only 12% [1]. Therefore the medical world nowadays aims towards alternative strategies such as thermotherapy so to increase this survival rate. Thermotherapy is a recently developed technique that consists in applying heat in the tumoral tissue. Radiofrequency ablation (RFA) heats the tissue through the use of one or multiple electrodes injecting radiofrequent currents. This technique is not yet well understood and engineering research is needed so to increase the efficacy of this method. This dissertation presents numerical techniques for the simulation of RFA. The configuration is bipolar in the sense that two needles are used for the injection of currents (Fig. 1). The mathematical model of RFA consists of three subproblems. First, a Poisson differential equation needs to be solved so to simulate the electric fields originating from the currents that are injected at a certain place. Secondly, the thermal behavior is characterized by the so-called Pennes bioheat equation for given heat sources and heat sinks. Finally, the temperature leads to an ablation of the tumoral tissues. This damage is described by an Arrhenius model. In addition, these three subproblems are strongly coupled to each other.

T

Fig. 1. Schematic representation of a bipolar configuration for RFA.

A validation of the different subproblems is provided in the thesis. In order to obtain a stable numerical solution of the problem in time and space, a stability analysis needs to be performed as well. Moreover, the values of the material properties are not exactly known. Consequently, their impact on the spatiotemporal distribution of the temperature and the corresponding extent of the ablated region needs to be investigated. Finally, an improvement of the existing bipolar radiofrequency ablation techniques is presented. II. M ETHODS The electromagnetic phenomena in RFA are modelled by expressing the continuity equation resulting in the following Poisson equation: ∇·[σ (~r, T (~r, t)) ∇φ (~r, t)] = I0 δ (~r − ~r1 )−I0 δ (~r − ~r2 ) (1)

Here, φ is the electric potential and σ is the electric conductivity, as a function of space (~r) and time (t). I0 is the volumetric current density that is originating from the electrodes at places ~r1 and ~r2 . Additionally, we use a Dirichlet type boundary condition ~ = −∇φ generates heat for eq. (1). The resulting electric field E 2 ~ (~r, t) . Eq. (1) is solved numerically by qE = σ (~r, T (~r, t)) E the finite difference method where we use a central approximation of the derivatives on a regular cubic grid. This results in a linear system of equations that is solved using the Preconditioned Biconjugate Gradient Stabilized method [2]. RFA introduces qE as an additional heat source next to the biological heat transfer model of Pennes: ∂T (~r, t) ρ (~r) c (~r) = ∇ · [k (~r, T (~r, t)) ∇T (~r, t)] ∂t −α (~r, t) ωρb cb (T (~r, t) − Tc ) + qE (2) with temperature T , mass density of the tissue ρ (of blood ρb ), specific heat capacity of the tissue c (of blood cb ), thermal conductivity k, capillary blood perfusion rate ω and body temperature Tc . We solve eq. (2) with the modified midpoint method and an additional Gragg-smoothing [3], which discretises this equation in the time domain. In order to obtain a stable and accurate solution, we performed a Von Neumann stability analysis. This analysis results in an upper limit for the timestep. Moreover, we validated the numerical solutions of eq. (1) and (2) for specific situations which could be solved analytically. In order to obtain a correct interpretation for thermotherapy, we employed an Arrhenius model for the thermal damage caused by RFA:   Zt ∆Ea − α (~r, t) = exp − Ae R·T (~r,t0 ) dt0  (3) 0

σ(T)

E

Thermal subproblem T

eq. (2) k(T)

Thermal damage eq. (3)

α

Fig. 2. Overview of the overall numerical model of RFA.

0.8

0.6 0.5 0.4 0.3 0.2 0

10

20

30

40

50

60

70

80

Time of current injection [s]

Fig. 4. Survival rate as function of the current injection time in the center of the simulation environment.

Fig. 5 illustrates the thermal damage in a cross section when using a certain number of pulses. In the beginning, the ablation zones are restricted to the regions close to the needle. Due to the influence of an optimized switch-off time, a temperature build-up arises in the intermediate non-ablated zone: the ablation zones can overlap with each other. This overlap increases when the number of pulses increases. Our concept could potentially help to treat tumors with larger dimensions (order of cm).

III. R ESULTS AND D ISCUSSION y [voxelnumbers]

50

50

start

40 30 20

10th pulse 1

40

0.9

30

0.8

20

0.7 10

10

0.6 10 50

20

30

x [voxelnumbers]

40

10

50 50

20th pulse

20

30

x [voxelnumbers]

40

50 0.5

50th pulse 0.4

y [voxelnumbers]

Using the developed numerical techniques, we investigated the influence of the distinct material parameters through numerical experiments on a tessellated geometry of 101 × 101 × 101 cubic voxels with a spatial discretisation of 1 mm. Fig. 3 shows for instance the temperature profile along the line that connects the two needles. We observe that the bipolar RFA technique is very local, i.e. the temperature rise due to the injection of currents is constrained to a small region of the order of millimeters. This figure also illustrates the influence of inhomogeneity (with tumor compared to overall healthy tissues) of material parameters on the temperature distribution.

= 0.5 s = 1.0 s = 1.5 s = 2.0 s = 2.5 s = 3.0 s

0.7

y [voxelnumbers]

φ

τf τf τf τf τf τf

40

y [voxelnumbers]

Electromagnetic subproblem eq. (1)

I0

1 0.9

Survival rate [-]

with the gasconstant R, frequency factor A and ∆Ea as the activation energy barrier. The survival rate α can be interpreted as the concentration of living cells at a certain time compared to its initial concentration. Traditionally, α in eq. (2) equals 1, whereas we modulate the Pennes bioheat equation using α since the effect of blood perfusion decreases when the tissue is thermally damaged. Fig. 2 illustrates how the above equations are entangled to each other resulting in a complex overall numerical model.

30 20 10

10

20

30

x [voxelnumbers]

40

50

40 0.3 30

0.2

20

0.1

10

0

10

20

30

x [voxelnumbers]

40

50

Fig. 5. Spatial distribution of the survival rate in the cross section z = 26.

IV. C ONCLUSION without tumor with tumor

350

Temperature [K]

345 340 335 330 325 320 315 310 42

44

46

48

50

52

54

56

58

60

y [voxelnumbers]

Fig. 3. Comparison of a temperature profile in the case of a geometry with a tumor and without a tumor.

Traditionally, the volumetric current density I0 is injected until the maximum temperature in the tissue reaches 80 ◦ C. The current is switched on again on an arbitrary moment. This leads to a local damage of the tumoral tissue. We propose the use of pulsed currents with fixed switch-off times. The longer the switch-off time, the wider the temperature distribution in space and thus the larger the ablation region. However, the absolute temperature will decrease as well resulting in less thermal damage. Consequently, a trade-off is present. Our numerical methods allow to calculate the influence of the different switch-off times on the thermal damage for fixed needle positions. Large differences are present between the various switch-off times as depicted in Fig. 4. Optimizing this switch-off time thus contributes to a more efficient treatment.

This work presents numerical methods for the interpretation of RFA. We use a finite difference method for the electromagnetic problem, whereas a timestepping method is used for the thermal problem. Next to a validation of the model, we also perform a stability analysis that provides a restriction on the timesteps. Based on the numerical results, bipolar RFA gives rise to a very local ablation region. Improvement of RFA is provided by using pulsed currents with fixed switch-off times. The developed methods are able to efficiently calculate the phenomena in RFA for given patient geometries, material parameters and needle positions. These methods can be used to increase the efficacy of RFA in the future. R EFERENCES [1] H. B. El-Serag, “Hepatocellular carcinoma,” New England Journal of Medicine, vol. 365, no. 12, pp. 1118 – 1127, 2011. [2] N. De Geeter, G. Crevecoeur, and L. Dupré, “An efficient 3-D eddy-current solver using an independent impedance method for transcranial magnetic stimulation,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 2, pp. 310 – 320, 2011. [3] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++: The Art of Scientific Computing, Cambridge University Press, 2002.

Inhoudsopgave Voorwoord

iv

Overzicht

vi

Extended Abstract

vii

Lijst van afkortingen

xii

1 Algemene inleiding 1.1 Motivatie . . . . . . . 1.2 Probleemstelling . . . 1.3 Doelstelling . . . . . . 1.4 Werkwijze en overzicht

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

2 Radiofrequente ablatie voor levertumoren 2.1 Thermotherapie . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Verschillende vormen van thermotherapie . . . . . . . 2.1.1.1 Radiofrequente ablatie . . . . . . . . . . . . . 2.1.1.2 Microgolfablatie . . . . . . . . . . . . . . . . 2.1.1.3 Cryoablatie . . . . . . . . . . . . . . . . . . . 2.1.2 Biologische warmteoverdracht . . . . . . . . . . . . . . 2.1.2.1 Pennes bio-warmtevergelijking . . . . . . . . 2.1.2.2 Weinbaum-Jiji bio-warmtevergelijking . . . . 2.1.2.3 De nieuwe aangepaste bio-warmtevergelijking 2.1.3 Biologische effecten van thermotherapie . . . . . . . . 2.2 Lever . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Anatomie . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Tumoren . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2.1 Metastasen . . . . . . . . . . . . . . . . . . . 2.2.2.2 Hepatocellulaire carcinoma . . . . . . . . . . 2.2.3 Behandelingsmethoden voor levertumoren . . . . . . . 2.2.4 Materiaaleigenschappen van levercellen . . . . . . . . 2.2.4.1 Elektrische geleidbaarheid . . . . . . . . . . . 2.2.4.2 Thermische geleidbaarheid . . . . . . . . . .

ix

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . . . . . . . .

. . . .

1 1 2 2 3

. . . . . . . . . . . . . . . . . . .

4 4 4 4 5 6 6 7 9 10 11 11 11 13 13 13 14 14 14 18

x

Inhoudsopgave 2.2.4.3

2.3

Overzicht van materiaaleigenschappen van ligne levercellen . . . . . . . . . . . . . . . Radiofrequente ablatie voor levertumoren . . . . . . . . . . 2.3.1 Actuele situatie . . . . . . . . . . . . . . . . . . . . . 2.3.1.1 Methode . . . . . . . . . . . . . . . . . . . 2.3.1.2 Complicaties . . . . . . . . . . . . . . . . . 2.3.1.3 Resultaten . . . . . . . . . . . . . . . . . . 2.3.2 Nieuwe configuratie . . . . . . . . . . . . . . . . . .

gezonde . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

en . . . . . . . . . . . . . .

ma. . . . . . . . . . . . . . . . . . . . .

19 20 20 20 21 21 22

3 Numerieke methoden voor radiofrequente ablatie 3.1 Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Discretisatie van simulatieomgeving . . . . . . . . . . . . . . . . . . . . . . . 3.3 Elektromagnetisch deelprobleem . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Wiskundige modellering . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Numerieke methoden voor het elektromagnetisch deelprobleem . . . . 3.4 Thermisch deelprobleem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Wiskundige modellering . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Modified midpoint-methode . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2.1 Beschrijving . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2.2 Numerieke fouten en stabiliteit . . . . . . . . . . . . . . . . . 3.5 Interpretatie van resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Thermische kostenfunctie . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1.1 Wiskundige modellering . . . . . . . . . . . . . . . . . . . . . 3.5.1.2 Numerieke methode . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Statistische verwerking . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Modulatie van de Pennes bio-warmtevergelijking . . . . . . . . . . . . . . . . 3.7 Overzicht van de globale werking van het numeriek model voor radiofrequente ablatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24 24 25 27 27 29 31 31 31 31 33 39 39 39 42 43 45

4 Numerieke analyse van numerieke methoden 4.1 Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Elektromagnetisch deelprobleem . . . . . . . . . . . . . 4.2.1 Analytische limietsituatie . . . . . . . . . . . . . 4.2.1.1 Delta-bron . . . . . . . . . . . . . . . . 4.2.1.2 Bolkapelektrode . . . . . . . . . . . . . 4.2.2 Numerieke simulatie . . . . . . . . . . . . . . . . 4.2.2.1 Overeenstemming met delta-bron . . . 4.2.2.2 Overeenstemming met bolkapelektrode 4.3 Thermisch deelprobleem . . . . . . . . . . . . . . . . . . 4.3.1 Analytische limietsituatie . . . . . . . . . . . . . 4.3.2 Numerieke simulatie . . . . . . . . . . . . . . . .

48 48 48 48 48 52 55 56 57 59 59 62

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

46

xi

Inhoudsopgave 5 Simulatie van thermotherapie 5.1 Configuratie van de simulatie . . . . . . . . . . . 5.2 Invloed van thermische en elektrische parameters 5.3 Invloed van discretisatie . . . . . . . . . . . . . . 5.4 Interpretatie van de numerieke simulatie . . . . .

. . . .

68 68 70 74 76

. . . . . . . . .

80 80 81 81 81 81 83 83 84 86

7 Conclusie 7.1 Overzicht van de masterproef . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Verder onderzoek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93 93 94

A De warmtevergelijking A.1 Afleiding van de warmtevergelijking . . . . . . . . . . . . . . . . . . . . . . . A.2 Fundamentele oplossing van de warmtevergelijking in een homogeen medium

96 96 97

. . . .

. . . .

6 Verbeterde thermotherapie voor leverkanker 6.1 Inleiding . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Gepulste thermotherapie . . . . . . . . . . . . . . . . 6.2.1 Simulatieomgeving . . . . . . . . . . . . . . . 6.2.2 Referentiestroomvorm . . . . . . . . . . . . . 6.2.2.1 Numerieke simulatie . . . . . . . . . 6.3 Verbeterde thermotherapie . . . . . . . . . . . . . . 6.3.1 Concept voor een efficiëntere thermotherapie 6.3.2 Aanschakeltemperatuur . . . . . . . . . . . . 6.3.3 Vaste afschakeltijd . . . . . . . . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

. . . .

. . . . . . . . .

Bibliografie

101

Lijst van tabellen

105

Lijst van figuren

106

Lijst van afkortingen BiCGSTAB Biconjugate Gradient Stabilized Method DNA desoxyribonucle¨ınezuur ECM extracellulaire matrix FN vals negatief FP vals positief FPF vals positieve fractie HCC hepatocellulaire carcinoma MRI Magnetic Resonance Imaging RFA radiofrequente ablatie TN waar negatief TP waar positief TPF waar positieve fractie

xii

Hoofdstuk 1

Algemene inleiding Cancer is a word, not a sentence. John Diamond

1.1

Motivatie

Kanker. Weinig woorden hebben dezelfde kracht om een onmiddellijke stilte af te dwingen. In 2008 werden wereldwijd 12.7 miljoen mensen geconfronteerd met de gevolgen van onbeheerste celdeling. 7.6 miljoen mensen stierven dat jaar aan de gevolgen van kanker. De kankerincidentie neemt nog voortdurend toe omdat de wereldpopulatie toeneemt, omdat de levensduur toeneemt en ook omdat kankerveroorzakend gedrag, zoals roken, toeneemt [28]. Conventionele therapieën van kanker zijn een combinatie van chirurgie, chemotherapie en radiotherapie. Helaas bieden ze niet bij alle vormen van kanker soelaas en zorgen ze vaak voor belastende nevenwerkingen die hun gebruik beperken. In deze thesis zullen we focussen op de behandeling van leverkanker. Leverkanker omhelst zowel de primaire levertumoren (hepatocellulaire carcinoma (HCC)), dit zijn kankers die hun oorsprong vinden in de levercellen zelf, als de secundaire levertumoren, de tumoren die elders in het lichaam ontstaan maar door metastasering de lever bereiken en zich daar verder ontwikkelen. Leverkanker is één van de meest voorkomende vormen van kanker. De incidentie van HCC bedraagt 740000 gevallen per jaar. Jaarlijks zijn er 694000 overlijdens door primaire leverkanker. De vijfjaarsoverleving bedraagt slechts 12% [20]. Chirurgie bij leverkanker is wel degelijk een effectieve behandelstrategie, maar kan slechts bij een beperkt aantal patiënten toegepast worden. Levertransplantaties zijn slechts van marginaal belang bij de behandeling van leverkanker gezien de beperkte beschikbaarheid van transplantlevers. Chemotherapie en radiotherapie blijken onvoldoende effectief te zijn. De medische wereld gebruikt thans in toenemende mate alternatieve strategieën, zoals thermotherapie, voor de behandeling van leverkanker en dit in een poging om de vijfjaarsoverleving te verhogen.

1

Hoofdstuk 1. Algemene inleiding

1.2

2

Probleemstelling

Thermotherapie bij kanker is een verzamelnaam voor behandelingsmethoden die de tumoren met warmte bestrijden. Zogenaamde globale thermotherapie warmt het volledige lichaam op en wordt toegepast wanneer een groot deel van het lichaam aangetast is door metastasen. In deze masterproef verrichten we onderzoek naar lokale thermotherapie van leverkanker onder de vorm van radiofrequente ablatie (RFA). Radiofrequente ablatie van leverkanker bestaat erin om het leverweefsel lokaal op te warmen door elektrische stromen te injecteren met behulp van één of meerdere elektroden. RFA is nog volop in ontwikkeling. Klinisch onderzoek heeft aangetoond dat de techniek suboptimaal is in de behandeling van leverkanker. Hier is een cruciale rol weggelegd voor ingenieurs. Om de techniek op punt te stellen dienen quasi-statische elektromagnetische veldberekeningen, die gekoppeld zijn met thermische vergelijkingen, uitgevoerd te worden. Op die manier is het mogelijk om een beter inzicht te verkrijgen in de lokale werking van RFA. Het is belangrijk dat een nauwkeurig numeriek model opgebouwd wordt dat de temperatuur in het weefsel kan berekenen, gegeven de geometrie van de lever en de tumor. Een levertumor kan immers ‘thermisch’ weggesneden worden indien de temperatuur ter hoogte van de tumor meer dan 50 ◦ C bedraagt. Temperaturen boven 100 ◦ C dienen vermeden te worden. Uiteraard zal de temperatuursdistributie afhankelijk zijn van de opbouw van de lever: de lever bestaat niet enkel uit levercellen maar ook uit andere weefsels. Bovendien bestaat er ook een grote onzekerheid omtrent de thermische en elektrische materiaaleigenschappen van deze biologische weefsels. Aangezien we nooit een exacte inschatting kunnen maken van elk weefsel of proces in de lever, moeten we onze toevlucht nemen tot wiskundige modellen. Het is primordiaal om een nauwkeurige inschatting te maken van de benaderingen inherent aan deze wiskundige modellen, teneinde een accuraat numeriek model op te stellen. Men is nu reeds in staat om tumoren met een diameter van 2 à 3 cm effectief te behandelen via een monopolaire naaldconfiguratie (één naald en één aardingsplaat). De grote uitdaging van thermotherapie is om grotere levertumoren efficiënt te behandelen. Het is noodzakelijk om numeriek na te gaan of dit mogelijk is via een bipolaire naaldconfiguratie (twee naalden) en wat de specificaties zijn voor de therapieparameters zoals de stroomvorm die door de naalden moet vloeien.

1.3

Doelstelling

De overkoepelende doelstelling van deze masterproef is om een numeriek model te ontwikkelen dat de elektromagnetische velden en de daaraan gekoppelde temperatuursprofielen simuleert. Het wiskundig model bestaat uit twee partiële differentiaalvergelijkingen met een sterke koppeling tussen beide vergelijkingen. Het oplossingsdomein is heterogeen, de materiaalparameters zijn temperatuursafhankelijk en de brontermen in het wiskundig model kunnen veranderen in de tijd en de ruimte. Voor een numeriek probleem van een dergelijke complexiteit dienen we er ons vooreerst van te vergewissen dat de oplossingsmethode stabiel is. Daartoe moet een grondige wiskundige analyse van de oplossingsmethode uitgevoerd worden.

Hoofdstuk 1. Algemene inleiding

3

Ten tweede dient nagegaan te worden wat de impact is van de materiaalparameters op de spatio-temporele temperatuursdistributie. Aangezien deze waarden onzekerheden zijn in het model, is het primordiaal om hun impact op de temperatuursdistributie te bepalen. Op die manier is het ook mogelijk om inzicht te verwerven in de werking van de techniek en tot hoe ver het lokale karakter van de techniek reikt. Ten derde wensen we na te gaan of de efficiëntie van de huidige techniek verbeterd kan worden.

1.4

Werkwijze en overzicht

Eerst wordt een literatuurstudie uitgevoerd in het kader van de geformuleerde probleem- en doelstellingen (hoofdstuk 2). Vertrekkende van die literatuurstudie worden in hoofdstuk 3 de verschillende deelproblemen van thermotherapie voor levertumoren gedefinieerd die zowel wiskundig als numeriek onderzocht dienen te worden. Om een overzichtelijke analyse van het probleem uit te kunnen voeren, opteren we voor een modulair softwareplatform. Op die manier is het eenvoudig om een welbepaalde geometrie (gezond en tumorweefsel, bloedvaten, etc.), naaldconfiguraties, stroomvormen, etc., in te geven. Dit laat ook toe dat dit platform na het proefschrift verder gebruikt en doorontwikkeld kan worden. De programmatie zelf gebeurt in een C++ omgeving. Verder zullen we in hoofdstuk 3 eveneens een kritische wiskundige analyse aanbrengen om de stabiliteit van de opgebouwde numerieke methoden te garanderen. Deze analyse laat toe om de benodigde discretisatieparameters te bepalen voor het bekomen van een stabiele oplossing van de wiskundige vraagstukken. Dit dient ge¨ımplementeerd te worden in het softwareplatform vooraleer de eigenlijke vergelijkingen opgelost worden. Op die manier is de numerieke methode universeel en dienen de discretisatieparameters niet door de gebruiker zelf bepaald te worden. Aangezien we RFA willen toepassen op menselijke levers, dienen we een extensieve validatie uit te voeren van ons numeriek model. In hoofdstuk 4 zullen we hiervoor de nauwkeurigheid van de verschillende deelproblemen apart analyseren. Voor de verwerking van de data resulterend uit de C++ omgeving, kunnen we Matlab gebruiken. Deze resultaten van onze numerieke methoden zullen daarna bediscussieerd worden in hoofdstuk 5. Enige aandacht gaat uit naar de invloed van de verschillende materiaalparameters teneinde inzicht te verwerven in de therapiemodaliteiten. Op basis hiervan geven we in hoofdstuk 6 mogelijke verbeteringen aan voor lokale thermotherapie. We besluiten deze masterproef in hoofdstuk 7 door de belangrijkste conclusies op te sommen, gevolgd door een vooruitblik naar mogelijk verder onderzoek.

Hoofdstuk 2

Radiofrequente ablatie voor levertumoren Mater artium necessitas.

2.1

Thermotherapie

Ter bestrijding van kanker, heeft men in de loop der jaren verschillende soorten behandelingsmethoden ontwikkeld. Recentelijk wordt ook thermotherapie naar voor geschoven als behandelingsmethode. Thermotherapie zal met behulp van warmte de kankercellen proberen te vernietigen. Naargelang de wijze waarop de warmte gegenereerd wordt in de buurt van de kankercellen onderscheiden we microgolfablatie, cryoablatie en radiofrequente ablatie (RFA).

2.1.1 2.1.1.1

Verschillende vormen van thermotherapie Radiofrequente ablatie

RFA maakt gebruik van een elektromagnetische golf met een frequentie in het radiofrequente gebied om een doelvolume te devitaliseren zonder het uit het lichaam te halen (virtuele chirurgische ablatie). RFA poogt met andere woorden biologische weefsels te vernietigen via elektriciteit van een ongemoduleerde sinuso¨ıdale golf met een elektromagnetische frequentie die karakteristiek is voor radiosignalen [35]. De meeste huidige RFA toestellen zijn monopolair, er is slechts één actieve elektrode en de elektrische stroom wordt gedissipeerd ter hoogte van de aarding. Voor de eigenlijke RFAprocedure wordt de elektrode in het doelvolume gebracht. Radiofrequente ablatie is dus een minimaal invasieve techniek: de incisie die nodig is om de elektrode ter hoogte van het tumorweefsel te brengen is veel kleiner vergeleken met deze van conventionele chirurgische ingrepen. Tijdens de RFA stroomt de alternerende radiofrequente stroom van de generator doorheen de geleidende tip van de elektrode en vervolgt zijn natuurlijk elektrisch pad tot zijn aarding. Zo ontstaat in feite een elektrisch circuit [35]. De ionen in het weefsel zullen trachten de sinuso¨ıdale verandering van de alternerende stroom te volgen. Dit heeft als gevolg dat er wrijvingswarmte ontstaat in het doelvolume (Ohmse verliezen) [11]. Ohmse verliezen worden 4

Hoofdstuk 2. Radiofrequente ablatie voor levertumoren

5

ook wel I 2 R-verliezen genoemd aangezien de warmte die ontstaat lineair evenredig is met de weerstand van het omgevende weefsel en kwadratisch evenredig met de lokale stroomdichtheid. Deze stroomdichtheid neemt kwadratisch af met de afstand tot de elektrode. Bijgevolg neemt de warmtegeneratie in het weefsel af met de afstand tot de elektrode tot de vierde macht. Aldus wordt het maximale vermogen afgezet in een heel klein gebied in de directe omgeving van de elektrode. De temperatuur in het doelvolume zal dus logischerwijze snel afnemen met toenemende afstand tot de elektrode [11, 52]. Het is de bedoeling van kankertherapie zoveel mogelijk maligne cellen te vernietigen en tezelfdertijd zoveel mogelijk gezonde cellen te sparen. Hier biedt RFA grote voordelen ten opzichte van de traditionele radiotherapie en chemotherapie. RFA biedt vooreerst de mogelijkheid om een heel lokale opwarming te veroorzaken in de buurt van de elektrode in tegenstelling tot bijvoorbeeld chemotherapie wat steeds een globale werking heeft op het lichaam. Verder heeft tumorweefsel ook andere fysische eigenschappen ten opzichte van normaal weefsel (zie subsectie 2.2.4). Hierdoor laat RFA toe om gebieden weefselspecifiek op te warmen. De opzet is om kankercellen een hogere temperatuur te laten bereiken dan de omliggende gezonde cellen en zo de tumor te vernietigen zonder al te veel schade aan te richten bij de gezonde cellen. De basis van de warmtegeneratie is zoals eerder vermeld de reactie van de ionen in het weefsel op de alternerende stroom. De frequentie van deze stroom dient radiofrequent te zijn. Radiogolven hebben een elektromagnetische frequentie gaande van 3 kHz tot 300 GHz. Meestal gebruikt men een frequentie in de buurt van 500 kHz. Deze frequentie is hoog genoeg om moleculaire wrijvingswarmte te veroorzaken zonder neuromusculaire reacties die optreden bij frequenties lager dan 20 kHz. De frequentie dient ook lager te zijn dan 20 MHz opdat de energietransmissie gebeurt zonder al te veel straling te veroorzaken [35]. De belangrijkste reden waarom een frequentie van 500 kHz wordt gekozen is omdat deze frequentie traditioneel gebruikt werd in bestaande elektrochirurgische toestellen [27]. Desalniettemin gaan er stemmen op om een lagere frequentie te hanteren [26, 25] (zie ook subsectie 2.2.4). 2.1.1.2

Microgolfablatie

Microgolfablatie is een andere vorm van thermotherapie waarbij de warmte ontstaat met behulp van microgolven. Deze hoogfrequente elektromagnetische golven worden gegenereerd door een antenne waardoor er een microgolf energieveld ontstaat binnen een eindig volume. Binnen dit energieveld zullen de watermoleculen in het weefsel roteren met de variërende elektrische velden waardoor opnieuw wrijvingswarmte ontstaat [36]. In tegenstelling tot de Ohmse verliezen bij RFA gaat het bij microgolfablatie dus eerder om diëlektrische verliezen [52]. Deze techniek ondervindt echter nog verschillende technologische problemen zoals bijvoorbeeld een impedantie mismatch tussen de antenne en het omgevende medium waardoor er een additionele warmtegeneratie ontstaat langs de voedingslijn en de efficiëntie van de antenne wordt aangetast [36]. De belangrijkste limitatie van microgolfablatie is echter dat de energie geabsorbeerd wordt in een zeer nauwe regio rond de antenne. De watermoleculen zullen de microgolven dermate sterk absorberen waardoor de penetratie van de elektromagnetische energie beperkt wordt tot 1 ` a 2 cm [52]. Desalniettemin verwacht men dat microgolfablatie


6

aan populariteit zal winnen naargelang de technologie aan maturiteit wint [36]. 2.1.1.3

Cryoablatie

Cryoablatie is een medische procedure waarbij ongewenst weefsel bevroren wordt aan de hand van een probe met een metalen tip, gekoeld door vloeibare stikstof. Hierna laat men het weefsel in situ ontdooien [3]. Celdood treedt op door het snel bevriezen van intracellulair water en de daaropvolgende lyse van de cel (het breken van het celmembraan) [36]. Bij het gebruik van warmte zoals in de vorige twee behandelingsmethoden blijken de maligne cellen gevoeliger te zijn aan letale schade dan gezonde cellen. Eigenaardig genoeg blijken tumorcellen meer resistent tegen letale schade door bevriezing in vergelijking met gezonde cellen [11]. Het is wel zo dat de vorm van tumoren significant verandert na cryoablatie. De disruptie van de tumorstructuur en de dehydratie zijn echter niet zo groot ten opzichte van normale levercellen [3]. Bijkomende complicaties van cryoablatie zijn onder andere het optreden van het cryoshock syndroom. Normaliter zal bloed enkel stollen indien er schade is aan de wand van het omhullende bloedvat. Bij het cryoshock syndroom ontstaat er echter een veralgemeende intravasculaire stolling. Zonder schade aan het bloedvat zal er toch een algemene stolling ontstaan van het bloed. Bijgevolg zullen meerdere organen falen. De oorzaak van dit syndroom is waarschijnlijk de vrijgave van cytokinines na het bevriezen. Deze complicatie heeft het gebruik van cryoablatie ernstig beperkt [36].

2.1.2

Biologische warmteoverdracht

Uit de vorige uiteenzetting is het duidelijk dat de diverse vormen van thermotherapie warmte injecteren ter hoogte van het tumorweefsel. Uit onze dagelijkse confrontatie met warmte weten we dat warmte niet gelokaliseerd blijft in de tijd. Warmte spreidt uit en geeft in verschillende media aanleiding tot verschillende temperatuursveranderingen. Bovendien wordt de warmteoverdracht in biologische weefsels complexer door de aanwezigheid van het thermoregulatoir systeem. Een studie van de biologische warmteoverdracht is dus essentieel om een modellering van thermotherapie mogelijk te maken. Een van de meest opmerkelijke eigenschappen van het menselijk thermoregulatoir systeem is zijn vermogen om de lichaamstemperatuur op circa 37 ◦ C te handhaven quasi onafhankelijk van de omgevingscondities. Dit wordt veroorzaakt door het enorme netwerk van bloedvaten in het menselijk lichaam. Bloed heeft een dubbele invloed op de thermische balans. Vooreerst is het in staat om als warmtebron of als koellichaam te ageren afhankelijk van de lokale weefseltemperatuur. Anderzijds laat de bloedstroom ook toe om de warmtedissipatie van het binnenste van het lichaam naar de omgeving te verhogen om tot een constante lichaamstemperatuur te komen [52]. Tot 1980 dacht men dat de warmteoverdracht gebeurde ter hoogte van de capillairen. Deze veronderstelling was gebaseerd op het feit dat enkel in de haarvaten zuurstof, CO2 , voedingsstoffen en afvalstoffen worden uitgewisseld. Bovendien heeft het haarvatennetwerk een zeer grote oppervlakte en dus veel mogelijkheden om warmte vrij te geven. Verschillende experimenten, uitgevoerd begin jaren 0 80, ontkrachten deze hypothese [7, 9].


7

Deze experimenten pogen de thermisch belangrijke bloedvaten te bepalen. Een belangrijke parameter in deze analyses is de thermische evenwichtslengte. De thermische evenwichtslengte van een bloedvat wordt gedefinieerd als de lengte waarover de temperatuur in het bloedvat afneemt met een factor e [9]. Indien deze thermische evenwichtslengte ongeveer dezelfde is als de fysische lengte van het bloedvat, wordt het bloedvat als thermisch belangrijk gecatalogeerd. Het resultaat van deze experimentele studies toont aan dat als de diameter van het bloedvat kleiner is dan 50 µm (d < 50 µm) de thermische evenwichtslengte veel kleiner is dan de fysische lengte waardoor alle warmtetransfer voor de intrede in deze vaten gebeurt. Indien d > 300 µm is de thermische evenwichtslengte veel groter dan de fysische lengte en treedt er dus quasi geen warmteverlies op. In het intermediaire gebied (50 µm < d < 300 µm) is de thermische evenwichtslengte ongeveer gelijk aan de bloedvatlengte waardoor we deze bloedvaten als thermisch significant kunnen beschouwen [7, 9, 52]. Ondanks de vooruitgang in het begrijpen van de biologische warmteoverdracht is het geen triviale opdracht om de warmteoverdracht in het menselijk lichaam te modelleren. De volgende drie redenen bemoeilijken de evaluatie van de warmteoverdracht [52]. • Het vasculair systeem is heel complex. Het is quasi onmogelijk om een numeriek model te ontwikkelen dat rekening houdt met alle thermisch belangrijke bloedvaten. • De thermische respons van bloedvaten hangt ook af van hun diameter en bloedperfusietempo. Deze hangen op hun beurt af van lokale parameters zoals de pH-waarde, O2 - en CO2 -concentraties en de temperatuur zelf. De bepaling van al deze externe en interne effecten op elk ogenblik en op elke positie is onmogelijk. • De thermische lengteschaal van thermisch belangrijke bloedvaten bedraagt zoals eerder vermeld ongeveer 50 − 300 µm. Deze temperatuurfluctuaties op micrometerschaal kunnen niet geresolveerd worden door de huidige technologie van toestellen die de temperatuur opmeten. Om de complexiteit van het vasculair systeem te omzeilen, opteert men in de literatuur vaak voor een continu¨ ummodel. Hier wordt de bloedstroom uitgemiddeld over een controlevolume en houdt men dus geen rekening met de eindige afmetingen van de bloedvaten. Het effect van de bloeddoorstroming wordt gemodelleerd door een additionele term in de warmtevergelijking van het weefsel [52]. Het verlies van detail omtrent het vasculair systeem resulteert uiteraard in het feit dat er geen uitspraken kunnen gedaan worden over de microscopische temperatuursvariaties. Meestal maakt men gebruik van het continu¨ ummodel afgeleid door Pennes [37]. 2.1.2.1

Pennes bio-warmtevergelijking

Bij de modellering van de biologische warmteoverdracht vertrekt Pennes van de continu¨ıteitsvergelijking die het behoud van energie uitdrukt voor een klein volume binnen het beschouwde lichaam. Zonder warmtebron resulteert de toepassing van behoud van energie binnen een klein volume in de gekende warmtevergelijking (zie sectie A.1). Elke bron van warmte zal een additionele term veroorzaken in de traditionele warmtevergelijking. Pennes beschouwde twee warmtebronnen in biologisch weefsel: enerzijds het effect van de bloedstroom die warmte aan-

8


of afvoert en anderzijds de metabolische warmte (de energie geproduceerd in elke cel) [37]. Wiskundig betekent dit dat: ∂Tt = ∇ · (kt ∇Tt ) + qbloed + qm (2.1) ∂t h i kg , c de specifieke warmtecapaciteit van Hierbij is ρ de massadichtheid van het weefsel m3 h i het weefsel kgJK , Tt de temperatuur van het weefsel ([K]), kt de thermische geleidbaarheid W en qm van het weefsel KWm , qbloed de term geassocieerd met bloed als warmtebron m 3 W de term geassocieerd met de metabolische warmteproductie m3 . ρc

afvoerader

venule

Quit

Qin

capillairen

Qin

Quit

toevoerslagader

arteriole

Figuur 2.1: Schematische voorstelling van een elementair volume weefsel dat slechts één slagaderader paar bevat.

Via fig. 2.1 kunnen we de sterkte van de bloedperfusieterm afleiden [52]. Hier wordt een schematisch diagram weergegeven van een klein volume weefsel dat één enkel slagader en ader paar bevat. De slagader staat via het haarvatennetwerk in contact met de ader. Indien men veronderstelt dat de temperatuur van de slagader en de ader constant blijft terwijl het bloed door dit weefselvolume stroomt, dan is de totale warmte die vrijgelaten is gelijk aan: Qρb cb (Ta − Tv ) = (Qin − Quit ) ρb cb (Ta − Tv )

(2.2)

Hierbij is Q de totale die het weefsel doorstroomt m3 h, ρb de i i massahbloedhoeveelheid kg J , cb de specifieke warmtecapaciteit van bloed , Ta de dichtheid van bloed kg K m3 temperatuur van de slagader ([K]) en Tv de temperatuur van de ader ([K]). Het volumetrisch warmteproductietempo qbloed is dan: qbloed =

(Qin − Quit ) ρb cb (Ta − Tv ) = ωρb cb (Ta − Tv ) V

(2.3)

ω stelt hierbij de hoeveelheid bloed voor die per volume-eenheid weefsel en per seconde door het weefsel stroomt s−1 [52]. In vgl. (2.3) is zowel Ta als Tv nog onbepaald. Pennes onderstelde in 1948 dat de biologische warmteoverdracht plaatsgreep ter hoogte van de capillairen


9

wegens hun grote oppervlakte. Aldus stelde hij de constant veronderstelde slagadertemperatuur Ta gelijk aan de lichaamstemperatuur Tc . Vervolgens benaderde hij de temperatuur van de ader Tv door de lokale weefseltemperatuur Tt aangezien het veneuze bloed in evenwicht treedt met het omringende weefsel ter hoogte van de haarvaten [52]. Met deze assumpties wordt qbloed uiteindelijk: qbloed = −ωρb cb (Tt − Tc ) (2.4) en bijgevolg krijgt de Pennes bio-warmtevergelijking de volgende vorm: ρc

∂Tt = ∇ · (kt ∇Tt ) − ωρb cb (Tt − Tc ) + qm ∂t

(2.5)

Een andere, meer wiskundige, manier om de basisveronderstellingen duidelijk te maken is door te stellen dat de temperatuur in de ader Tv een lineaire combinatie is van de slagadertemperatuur Ta en de weefseltemperatuur Tt [37]: Tv = Tt + κ (Ta − Tt )

(2.6)

Hierbij is κ de evenwichtsconstante die de afwijking ten opzichte van de weefseltemperatuur voorstelt (0 ≤ κ ≤ 1). Het aandeel van de bloedperfusie in de warmtevergelijking wordt dan: qbloed = ωρb cb (κ − 1) (Tt − Tc )

(2.7)

In de modellering van Pennes wordt κ = 0 aangezien Pennes veronderstelt dat Tv = Tt [37, 34]. Het grote voordeel van vgl. (2.5) is dat deze vergelijking relatief eenvoudig te gebruiken is en dat het de manipulatie van de bloedterm via twee parameters (ω en Tc ) toelaat. Bovendien zijn er behoorlijke overeenkomsten tussen theorie en experiment [37, 52]. De belangrijkste limitaties van dit model zijn verbonden met de basisassumpties. Zoals eerder vermeld grijpt de warmteoverdracht niet plaats ter hoogte van de capillairen [34]. De basisveronderstellingen dat Ta = Tc en Tv = Tt zijn bijgevolg foutief, immers: • De veronderstelling dat Ta constant blijft bij transport van het hart naar het bed van capillairen is enkel correct als de diameter van het bloedvat groter is dan 300 µm (wegens grote thermische evenwichtslengte). • Tv ≈ Tt is enkel correct als d < 50 µm (wegens kleine thermische evenwichtslengte is alle warmte al vrijgegeven vooraleer het bloed de ader bereikt) Men kan dus concluderen dat de Pennes bio-warmtevergelijking het effect van de bloedperfusie overschat [52]. Om tot een accurater model te komen moeten we rekening houden met enerzijds de temperatuursvariatie langs de slagader en anderzijds met de warmte die door het bloed in de aders, dat een tegengestelde stroomrichting heeft, terug wordt opgenomen [52]. 2.1.2.2

Weinbaum-Jiji bio-warmtevergelijking

Na experimenten begin de jaren 0 80, twijfelen meer en meer mensen aan de geldigheid van de Pennes bio-warmtevergelijking [7, 9]. Weinbaum en Jiji ontwikkelden een nieuwe bio-warmtevergelijking op basis van een anatomische analyse van een weefsel waarbij de dominante mode


10

van warmtetransfer de tegenstroom warmteuitwisseling was tussen een thermisch belangrijk slagader en aderpaar [48, 52]. Het bijna perfecte tegenstroom warmteuitwisselingsmechanisme impliceert dat de meeste warmte die de slagader verlaat opnieuw wordt opgenomen door de ader waarin het bloed in een tegenovergestelde richting stroomt. De warmte wordt dus niet vrijgegeven aan het weefsel in de directe omgeving. Ter vereenvoudiging van het probleem onderstelden Weinbaum en Jiji dat de tegenstroom warmteuitwisseling perfect was. Het resultaat is dat Ta = Tv of in termen van de evenwichtsconstante: κ = 1 [34]. De additionele bloedperfusieterm van vgl. (2.7) verdwijnt met andere woorden. Indien er een gradiënt bestaat van de weefseltemperatuur langsheen de assen van de ader en de slagader (∇Tt 6= 0), dan zullen de ader en de slagader verschillende hoeveelheden energie transfereren in een vlak loodrecht op hun assen zelfs al is er geen netto warmteuitwisseling in de capillairen [52]. Dit resulteert in een netto energie transfer die equivalent is aan een verhoging van de thermische geleidbaarheid van het weefsel: kef f = kt [1 + f (ω)] (2.8) Zoals aangeduid in vgl. (2.8) is deze effectieve thermische geleidbaarheid een functie van de bloedperfusiesnelheid ω. De bio-warmtevergelijking die resulteert uit de perfecte tegenstroom warmteuitwisseling wordt de Weinbaum-Jiji bio-warmtevergelijking genoemd: ∂Tt = ∇ · (kef f ∇Tt ) + qm (2.9) ∂t De beperking van het model volgens Weinbaum en Jiji is het feit dat de tegenstroom warmteuitwisseling maar dominant is in een beperkt aantal weefsels (voornamelijk perifere weefsels). In relatief grote bloedvaten (d > 200 µm) is het niet zo dat de meeste warmte die de slagader verlaat opnieuw wordt opgenomen door de nabije ader [48]. De Weinbaum-Jiji biowarmtevergelijking is dan ook het meest geschikt voor de beschrijving van warmteoverdracht in dichtgelegen paren van microbloedvaten in spierweefsel [52]. Bovendien is de bepaling van kef f complex en heeft men nood aan accurate anatomische gegevens die niet voor elk belangrijk bloedvatenpaar beschikbaar zijn. Dit in tegenstelling tot de vergelijking van Pennes waarbij enkel ω en Tc nodig zijn om het effect van de bloedperfusie te bepalen [52]. ρc

2.1.2.3

De nieuwe aangepaste bio-warmtevergelijking

Het is duidelijk dat het Pennes en het Weinbaum-Jiji model twee extreme situaties voorstellen wat betreft biologische warmteoverdracht. Waar het model volgens Pennes net veronderstelt dat het arterieel bloed alle warmte vrijgeeft aan het omgevend weefsel, onderstelt het Weinbaum-Jiji model dat alle warmte die de slagader verlaat opnieuw opgenomen wordt door de ader waarin het bloed in een tegenovergestelde richting stroomt. In termen van de evenwichtsconstante κ komt de extremiteit zo mogelijk nog meer tot uiting: κ = 0 in het Pennes model en κ = 1 in het Weinbaum-Jiji model. Enkele theoretici [49] stellen voor om het hoofd te bieden aan de tekortkomingen van beide modellen door een correctiefactor te hanteren in de bloedperfusieterm van het Pennes model: ∂Tt = ∇ · (kt ∇Tt ) − ωρb cb (Tt − Ta0 ) + qm (2.10) ∂t Hierbij is = 1 − κ. Deze correctiefactor zorgt ervoor dat de overschatting van de bloedperfusieterm in de originele Pennes bio-warmtevergelijking gecorrigeerd wordt door rekening ρc


11

te houden met het feit dat een zeker percentage van de warmte die door de slagaders wordt vrijgegeven opnieuw opgenomen wordt door de aders [52]. wordt op basis van empirische resultaten bepaald [34]. Het blijkt dat vooral afhangt van de afstand tussen de slagader en ader paren en de diameter van de bloedvaten [52]. Voor verschillende skeletspieren bedraagt de experimenteel bepaalde tussen 0.6 en 0.8 [52]. Dit toont aan dat een intermediaire toestand de biologische warmteoverdracht beter beschrijft dan de extreme modellen van Pennes en Weinbaum-Jiji. Helaas zijn voor de levercellen geen gegevens omtrent beschikbaar. Daarnaast gebruikt men in vgl. (2.10) vaak Ta0 in plaats van Tc . De veronderstelling dat Ta0 = Tc blijkt enkel te kloppen indien de bloedperfusiesnelheid vrij hoog is [52].

2.1.3

Biologische effecten van thermotherapie

De biologische effecten van de verschillende vormen van thermotherapie zoals RFA betreffen verschillende en complexe mechanismen. Deze mechanismen hangen enerzijds af van de temperatuur en de duur van de warmteblootstelling en anderzijds van lokale factoren zoals de orgaanperfusie, weefseldichtheid en elektrolytconcentraties. Naargelang de hoogte van de temperatuur zie je andere celschade optreden. Vanaf 42 ◦ C begint de thermische schade. Deze gaat gepaard met een progressieve degeneratie van de cellen door schade aan het celmembraan, de chromosomen (desoxyribonucle¨ınezuur (DNA)) en het cytoskelet. Daarnaast ontstaat er ook een vertraging van de metabolische werking van de cel en een verstoorde eiwitsynthese [35]. Naargelang de temperatuur stijgt zal de blootstellingstijd die nodig is om een irreversibele celschade aan te richten exponentieel afnemen. Vanaf een temperatuur van 60 ◦ C worden de intracellulaire prote¨ınen volledig gedenatureerd en vindt ook de vernietiging plaats van het celmembraan door het oplossen en smelten van de vetdubbellaag waaruit het bestaat [11]. Uiteraard zal een temperatuur boven de 100 ◦ C naast de eerdere fenomenen ook het koken en verdampen van intracellulair water met zich meebrengen. Dit resulteert in een quasi onmiddellijke celdood en men spreekt van coagulation necrosis [11]. We kunnen besluiten dat een temperatuur tussen 50 ◦ C en 100 ◦ C doorheen het targetvolume als doeltreffend wordt aanzien bij thermotherapie [35]. Naast celdood heeft de warmte ook een effect op de vasculatuur in het doelvolume. Na de blootstelling zullen verschillende histopathologische veranderingen ontstaan zoals intravasculaire tromboses en neutrofiele aanhechting aan het endotheel van de bloedvaten. Tijdens RFA neemt men inderdaad een afname waar in de perfusie door de microvasculatuur en na de RFA ontstaat er een vasculaire shutdown, i.e. de beschadigde bloedvaten laten geen doorstroming meer toe [35].

2.2 2.2.1

Lever Anatomie

De lever is het grootste orgaan in de buikholte. Ze weegt ongeveer 1.5 kg en ligt onder het middenrif, rechtsboven in de buikholte. De lever heeft twee grote intredende bloedvaten: het meeste bloed wordt aangevoerd door de poortader die van de darmvlokken komt (zuurstofarm bloed dat echter rijk is aan voedings-


12

stoffen) en de leverslagader die van de aorta komt (en zuurstofrijk bloed aanbrengt). Er zijn twee grote uitgaande bloedvaten: de linker- en de rechter leverader die het bloed komende van de lever afvoeren naar de onderste holle ader. Al het bloed dat uit de maag en de darm komt moet dus eerst de lever passeren. De bloedstroomsterkte door de lever is afhankelijk van het verschil in druk in de aanvoerende en afvoerende vaten en van de weerstand in het leverstroombed. Er is echter een dubbele aanvoer: lage druk in poortadersysteem (1 kPa) en hoge druk in leverslagader (12 kPa). De leverdoorbloeding is het sterkst na een maaltijd. Voor een volwassen man is dat 1.5 l per minuut na een maaltijd. Inwendig bestaat de lever uit duizenden leverlobjes die als veelhoekige kolommetjes van 1 ` a 2 mm naast elkaar liggen (fig. 2.2). Ieder lobje bevat centraal een vertakking van één van beide leveraders. Verder zijn er strengen van levercellen die radiaal gerangschikt zijn. Ze worden bevloeid met bloed dat de spleetvormige ruimtes tussen de strengen vult. Dat bloed is afkomstig van de poortader en van de leverslagader en vloeit af naar een vertakking van de leverader. In de wanden van de bloedruimten liggen de stercellen van Kupffer. De cellen van Kupffer fagocyteren de afstervende rode bloedcellen. Daarbij breken ze hemoglobine af tot ijzer en een geel pigment, bilirubine, die beide worden opgeslagen door de levercellen. Het ijzer wordt als reserve opgeslagen. De levercellen nemen uit het bloed dat erlangs stroomt stoffen op en geven er stofwisselingsproducten aan af.

Figuur 2.2: Schematische opbouw van een leverlobje [1].

Alleen de gal volgt een andere weg. Galsap wordt verzameld in spleetjes die tussen de levercellen liggen, de galkanaaltjes. Die kanaaltjes komen samen in de galgangen. Gal en bloed stromen in tegengestelde zin. De lever speelt een grote rol bij de homeostase, het constant houden van het inwendige milieu in het lichaam:


13

• De lever regelt de concentraties van de opgenomen voedingsstoffen en van verschillende bloedsubstanties zoals plasma-eiwitten en stollingsfactoren. • Ook bij de intermediaire stofwisseling speelt de lever een sleutelrol. Uit cholesterol worden galzouten aangemaakt die een rol spelen bij de vetvertering en de vetopname. • De lever heeft een belangrijke detoxificatiefunctie (voor zowel afvalstoffen uit de stofwisseling als lichaamsvreemde stoffen). • De lever is een opslagplaats voor vitaminen. Tot slot moet opgemerkt worden dat het zuurstofverbruik van de lever ongeveer 20% is van het totale lichaamsverbruik. Hierdoor levert de lever een grote rol bij de productie van lichaamswarmte (1700 kJ per dag).

2.2.2 2.2.2.1

Tumoren Metastasen

Metastasen ontstaan wanneer een kwaadaardige primaire tumor groeit tot in een bloedvat en/of een lymfevat. Cellen van deze kwaadaardige tumor komen los en worden meegevoerd met het bloed of de lymfe. De losgekomen tumorcellen gaan zich elders in het lichaam vastzetten en hier ontstaat een secundaire tumor of metastase. Metastasen komen vaak voor in de lever omdat de lever een bloedvatrijk orgaan is. Primaire tumoren die vaak metastaseren in de lever zijn: colorectale tumoren, maagtumoren, borstkanker, pancreascarcinoom, longcarcinoom. Levermetastasen geven aanvankelijk geen klinische symptomen. Wanneer er klachten zijn dan zijn ze vaak atypisch: opgezette lever, misselijkheid en braken, verminderde eetlust, vermoeidheid. Jeuk en geelzucht ontstaan wanneer de metastase de galwegen blokkeert. Grote levermetastasen leiden naar leverinsufficiëntie en uiteindelijk leverfalen: er ontstaat een hepatisch coma gevolgd door overlijden. 2.2.2.2

Hepatocellulaire carcinoma

Dit is de vijfde meest voorkomende kanker bij mannen en de zevende bij vrouwen [20]. Deze HCC komen vooral voor in ontwikkelingslanden. Mannen worden bovendien meer getroffen dan vrouwen. In 2008 waren er wereldwijd 694000 overlijdens tengevolge van HCC. HCC heeft een unieke geografische, geslachts- en leeftijdsverdeling die bepaald worden door zijn specifieke etiologische factoren. Omwille van zijn hoge fataliteit (verhouding mortaliteit ten opzichte van incidentie = 0.93) is leverkanker de derde meest voorkomende oorzaak van dood door kanker wereldwijd. HCC zijn agressieve tumoren die meestal ontstaan in het verloop van een chronische leverziekte of cirrose. Leverfibrose (voorstadium van levercirrose) is een vorm van wondheling waarbij beschadigde (door bijvoorbeeld een chronische virale infectie) levergebieden worden ingekapseld door een extracellulaire matrix (ECM) of littekenweefsel. De snelheid van evolutie naar fibrose hangt af van enerzijds de oorzaak van het leverlijden en anderzijds van gastfactoren.


14

Leverfibrose ontwikkelt zich over een periode van jaren. Aanvankelijk is deze evolutie omkeerbaar, progressieve fibrose leidt onherroepelijk naar cirrose. Momenteel hebben we nog geen antifibrotische geneesmiddelen ter onze beschikking. Bij een fibrotische lever zijn er significante veranderingen in de ECM: het collageengehalte neemt bijvoorbeeld toe met een factor 3 tot 10. Dit heeft een grote invloed op de functie van de levercellen en verklaart de synthetische en metabole dysfunctie geobserveerd in patiënten met gevorderde fibrose. De bron voor ECM zijn de leverstellaatcellen. Fibrose is een uiting van een onevenwicht in de balans tussen matrixproductie en afbraak. De diagnose wordt vaak laat gesteld in het verloop van de ziekte en de gemiddelde overleving na het stellen van de diagnose is beperkt (6 tot 20 maanden). De therapeutische opties worden bepaald door de hepatische reserve [45]. De belangrijke risicofactoren voor de ontwikkeling van HCC zijn hepatitis B-virus dragerschap, chronische hepatitis C-virus infectie, erfelijke hemochromatose en cirrose van om het even welke oorsprong [13].

2.2.3

Behandelingsmethoden voor levertumoren

De selectie van de therapie wordt bepaald door de ernst van het onderliggend leverlijden, de grootte en de distributie van de intrahepatische tumoren, de bloedtoevoer en de algemene toestand van de patiënt. Hieronder wordt een overzicht gegeven van de verschillende therapieën. Er zijn enerzijds chirurgische therapieën: • Resectie: vooral voor patiënten zonder cirrose met HCC in de beginstadia (5% van de patiënten) [20]. Hierbij is er trouwens een kans van 70% dat er nieuwe HCC verschijnen na 5 jaar [20]. • Transplantatie: deze mogelijkheid is vaak hypothetisch gelet op het beperkt aantal beschikbare transplantlevers. Anderzijds zijn er ook niet-chirurgische therapieën: • Systemisch: chemotherapie. Grootste nadeel is dat systemische behandelmethodes een globale werking hebben op het hele lichaam. • Lokaal: percutane alcoholinjectie, tie en RFA [20].

90 Y

-radioembolisatie, transarteriële chemoembolisa-

Ondanks de vele therapieën is men er vooralsnog niet in geslaagd om de hoge fataliteit te reduceren.

2.2.4 2.2.4.1

Materiaaleigenschappen van levercellen Elektrische geleidbaarheid

Frequentie-afhankelijkheid De frequentie-afhankelijkheid van de elektrische geleidbaarheid σ van een biologisch weefsel is belangrijk bij toepassing van RFA. Afhankelijk van de waarde van de elektrische geleidbaarheid zal het weefsel meer of minder opwarmen en dus


15

meer of minder schade ondervinden. De elektrische geleidbaarheid staat in verband met de complexe relatieve permittiviteit ˆ [22, 38]: ˆ = 0 − j00 σ 00 = ω0

(2.11)

ˆ is afhankelijk van de frequentie van het aangelegde elektrische veld ω. Er is immers altijd een vertraging tussen veranderingen in de (moleculaire) polarisatie van het diëlektricum en de verandering van het elektrisch veld. De eerste die een vergelijking opstelde voor de diëlektrische relaxatie was Debye [24, 38]: ˆ (ω) = ∞ +

s − ∞ 1 + jωτ

(2.12)

In de Debye-vergelijking stelt ∞ de relatieve permittiviteit voor gemeten bij hoge frequenties (ωτ >> 1) en s de relatieve permittiviteit gemeten bij lage frequenties (ωτ << 1). Verder is τ de karakteristieke relaxatietijd van het diëlektricum. Bijgevolg is ook 00 en dus σ frequentie-afhankelijk: 00 (ω) = =

ωτ (s − ∞ ) 1 + ω2τ 2 σ (ω) ω0

(2.13)

σ (ω) bevat meestal nog een frequentie-onafhankelijke component die overeenkomt met de statische ionische conductiviteit σ0 : σ (ω) = σ0 + σd (ω) [22, 24, 38]. Experimenten tonen aan dat de frequentie-afhankelijkheid van biologische weefsels niet helemaal verklaard wordt door de Debye-vergelijking (vgl. (2.12)) [24]. Biologische weefsels zijn dermate complex qua structuur en samenstelling dat elke dispersieregio kan verbreden door veelvuldige contributies in deze regio. Deze verbreding (distributie van relaxatietijden) kan empirisch in rekening gebracht worden door een distributieparameter α te gebruiken in de Debye-vergelijking. Aldus bekomen we de Cole-Cole vergelijking [24, 38, 10]: ˆ (ω) = ∞ +

s − ∞ 1 + (jωτ )1−α

(2.14)

Bovendien zal elk type polarisatie zijn eigen karakteristieke respons (eigen relaxatietijd) hebben tegenover het aangelegde elektrische veld [38]. Bij biologische weefsels zijn er drie dispersieregio’s te onderscheiden met elk een eigen relaxatietijd, nl. α-, β- en γ-dispersie [22, 38]. Deze dispersieregio’s worden weergegeven in fig. 2.3. De belangrijkste eigenschappen van het diëlektrisch spectrum van biologische weefsels, zoals weergegeven in fig. 2.3, zijn [22]: • De relatieve permittiviteit 0 neemt af met toenemende frequentie, terwijl de elektrische geleidbaarheid σ toeneemt met de frequentie. • Een eerste regio (α-dispersie) in het frequentiegebied van 10 Hz tot 1 kHz wordt geassocieerd met diffusieprocessen van ionen ter hoogte van het celmembraan. • Vervolgens is er β-dispersie tussen 100 kHz en 10 MHz die de polarisatie van het celmembraan zelf voorstelt. Andere contributies zijn de polarisatie van prote¨ınen en verschillende organische macromoleculen.

16


8

10

²0 σ

α 6

10

²0 [-] σ [S/m]

4

10

β

γ

2

10

0

10

−2

10

0

10

2

10

4

6

10

10

8

10

10

10

12

10

Frequentie [Hz]

Figuur 2.3: Diëlektrisch spectrum van biologische weefsels.

• Tenslotte is er een derde regio (γ-dispersie) voor frequenties groter dan 1 GHz. Deze regio wordt veroorzaakt door de polarisatie van watermoleculen. Het diëlektrisch spectrum van biologische weefsels kan met andere woorden beter beschreven worden door veelvuldige Cole-Cole dispersie: ˆ (ω) = ∞ +

X n

∆n 1 + (jωτn )1−αn

(2.15)

Hierbij is n het aantal dispersiemodes, bij biologisch weefsel is n = 3 (α, β, γ). Toch gebruikt men vaak het vier Cole-Cole model (n = 4) om meer flexibiliteit te verkrijgen bij het fitten van de theoretische vergelijking aan de experimentele data [24]. In de literatuur zijn heel wat gegevens beschikbaar over de permittiviteit en de elektrische geleidbaarheid van verschillende biologische weefsels waaronder de lever [22, 23, 24]. Wel is het zo dat er een grote variatie in de experimentele waarden zit omdat de meetmethode ex vivo of in vivo kan gebeuren en er niet altijd op menselijk weefsel wordt gemeten. Bij toepassing van RFA werken we in het frequentiegebied van 20 kHz tot 20 MHz (zie subsectie 2.1.1.1). Dit gebied komt ruwweg overeen met het gebied van β-dispersie. Het effect van β-dispersie wordt vooral veroorzaakt door de specifieke diëlektrische eigenschappen van de celmembranen. Celmembranen zijn samengesteld uit een vetdubbellaag die als een capaciteit kan gemodelleerd worden tussen de intra- en extra-cellulaire ruimte [26]. Bij frequenties beneden de β-relaxatiefrequentie wordt de stroom vooral beperkt tot de extra-cellulaire ruimte. Boven de β-relaxatiefrequentie laat de capacitieve koppeling van de celmembranen de stroom toe om door het membraan te gaan waardoor de elektrische geleidbaarheid van het weefsel stijgt [26]. Dit is onmiddellijk een intu¨ıtieve verklaring voor het feit dat de elektrische geleidbaarheid toeneemt met de frequentie. Meestal stelt men de frequentie in op de traditionele waarde, nl. 500 kHz [27]. Verschillende artikels in de literatuur stellen echter voor om een lagere frequentie te hanteren [26, 25]. Aangezien de eigenschappen van tumorcellen verschillend zijn van die van gezonde cellen


17

verwachten we dat de elektrische geleidbaarheid ook (lichtjes) verschillend zal zijn. D. Haemmerich et al. verzamelden in vivo resultaten van de elektrische conductiviteit bij 17 ratten (24 levertumoren) [26, 27, 25]. De resultaten van deze metingen worden weergegeven in fig. 2.4.

Figuur 2.4: Elektrische geleidbaarheid van normale levercellen en tumorcellen in functie van de frequentie [26].

Het is vooreerst duidelijk dat tumorweefsel een hogere geleidbaarheid heeft dan normale levercellen. Het verschil in elektrische geleidbaarheid wordt echter steeds kleiner naarmate de frequentie toeneemt. Ex vivo-studies van menselijke tumoren tonen daarenboven aan dat dit verschil eigenlijk nog groter is dan de in vivo-studie van rattumoren aantoont [25]. Voor therapeutische doeleinden is het interessant dat er een groot verschil is tussen de geleidbaarheden van maligne en gezonde cellen. Het temperatuursverschil tussen de tumorcellen en de gezonde cellen zal dan immers ook groter zijn, waardoor er meer schade kan aangericht worden aan de kwaadaardige cellen, terwijl de gezonde cellen relatief minder thermische schade ondervinden. Dit wordt bevestigd door een numerieke simulatie [25]. We kunnen dus concluderen dat voor de toepassing van RFA best een zo laag mogelijke frequentie gebruikt wordt. Deze frequentie dient echter groter te zijn dan 20 kHz om neuromusculaire reacties te vermijden [35]. Tot op heden zijn er echter nog geen studies bekend van radiofrequente ablatie bij een lage frequentie. Temperatuursafhankelijkheid Zoals eerder vermeld werkt RFA vooral in het gebied dat door β-dispersie gekenmerkt wordt. Beneden de β-relaxatiefrequentie wordt de stroom beperkt tot de extra-cellulaire ruimte. Bij toename van de temperatuur zal de viscositeit van de extra-cellulaire vloeistof echter afnemen waardoor de mobiliteit van de ionen die de stroom voeren toeneemt [33]. Dit heeft uiteraard als gevolg dat de elektrische geleidbaarheid zal toenemen met de temperatuur. Boven de β-relaxatiefrequentie zal een deel van de stroom in het extra-cellulaire gebied door het celmembraan diffunderen. Aangezien de diffusie toeneemt met toenemende temperatuur zal ook boven deze frequentie de elektrische geleidbaarheid toenemen. De elektrische geleidbaarheid blijkt zelfs lineair toe te nemen met de temperatuur


18

[33, 27]: σ (T ) = σ0 · [1 + σ1 (T − T0 )]

(2.16)

Afhankelijk van de bron bedraagt σ1 1.5% K−1 ([27]) tot 2% K−1 ([33]). Deze σ1 blijkt dezelfde te zijn voor tumorweefsel als voor gezonde cellen [27]. In dit model verwaarlozen we niet-lineaire fenomenen door irreversibele schade aan de weefsels [27]. Tot op heden zijn er geen gegevens beschikbaar om deze hypothese te rechtvaardigen of te ontkrachten. Intu¨ıtief voelen we aan dat deze veronderstelling gerechtvaardigd zal zijn als de maximumtemperatuur in het weefsel beperkt blijft (Tmax < 100 ◦ C). 2.2.4.2

Thermische geleidbaarheid

Temperatuursafhankelijkheid De temperatuursafhankelijkheid van de thermische geleidbaarheid k van leverweefsel is niet eenvoudig te bepalen. De eerste in vitro metingen geven aan dat de thermische geleidbaarheid van menselijk leverweefsel de volgende lineaire uitdrukking aanneemt [46]: k (T ) = 0.4692 + 0.001161 T (2.17) W Deze relatie tussen k Km en T ([◦ C]) is geldig in het temperatuursgebied tussen 3 ◦ C en 45 ◦ C. We kunnen dit ook schrijven als: k (T ) = 0.5122 + 0.001161 (T − 37 ◦ C)

(2.18)

waarbij k (T ) wordt uitgezet in functie van het relatief temperatuursverschil ten opzichte van de evenwichtstemperatuur van het biologisch weefsel (37 ◦ C). Biologische weefsels vertonen blijkbaar bijna een zelfde temperatuursafhankelijkheid als water waarbij de lineaire voorfactor 0.001575 bedraagt. Dit hoeft ons niet te verbazen aangezien water een belangrijk bestanddeel van biologisch weefsel is. Een groot probleem gepaard met in vitro experimenten is het feit dat biologische weefsels zullen denatureren als ze geen voedingsstoffen meer aangereikt krijgen. De samples worden dus op een onbepaalde manier veranderd ten opzichte van hun toestand in vivo [46]. Tot heden zijn er enkel in vivo metingen uitgevoerd op varkenslevers [2]. Deze resultaten geven weer dat de thermische geleidbaarheid in levend biologisch weefsel sterk toeneemt als gevolg van de bloedperfusie. Wij zijn echter meer ge¨ınteresseerd in de thermische geleidbaarheid van de levercellen an sich, aangezien het effect van de bloedperfusie door een aparte term wordt gemodelleerd in de Pennes bio-warmtevergelijking. Opnieuw verwaarlozen we niet-lineaire effecten door onomkeerbare schade aan de weefsels. In de literatuur is een experiment terug te vinden waarbij in vitro metingen worden uitgevoerd op de lever van een koe [2]. Hierbij laat men de temperatuur tot een bepaalde waarde toenemen en vervolgens koelt men het weefsel geleidelijk af tot zijn initiële temperatuur. Tot een maximumtemperatuur van 80 ◦ C komt de thermische geleidbaarheid van het opgewarmde en vervolgens afgekoelde sample overeen met de initiële thermische geleidbaarheid. Dit duidt aan dat er nog geen irreversibele schade is aangericht aan de levercellen en we de lineaire relatie tussen k en T mogen hanteren. Het temperatuursgebied waarin de lineaire relatie van vgl. (2.18) geldt mag met andere woorden uitgebreid worden tot 80 ◦ C. Bij een maximumtemperatuur van 90 ◦ C is dit niet langer het geval en is de thermische geleidbaarheid na de

19


volledige cyclus significant verschillend van deze voor de cyclus. Dit wordt weergegeven in fig. 2.5.

Figuur 2.5: Experimentele verwarmings- en afkoelingscyclus tot 90 ◦ C waarbij de relatie tussen k en T niet-lineaire effecten vertoont [2].

2.2.4.3

Overzicht van materiaaleigenschappen van gezonde en maligne levercellen

In het vervolg van deze masterproef gebruiken we RFA bij een frequentie van 500 kHz. De elektrische en thermische geleidbaarheden voor levercellen in het temperatuursgebied tot 80 ◦ C worden gegeven door lineaire relaties: σ (T ) = σ0 · [1 + σ1 (T − T0 )]

(2.19)

k (T ) = k0 + k1 (T − T0 )

(2.20)

De verschillende coëfficiënten voor gezonde en maligne levercellen vatten we samen in tabel 2.1 [27]. Tabel 2.1: Materiaaleigenschappen van levercellen bij toepassing van RFA [27].

Type weefsel Gezond leverweefsel Maligne leverweefsel

σ0 S m−1 0.36 0.45

σ1 K−1 0.015 0.015

k0 W m−1 K−1 0.53 0.53

k1 W m−1 K−2 0.00116 0.00116

20


2.3 2.3.1 2.3.1.1

Radiofrequente ablatie voor levertumoren Actuele situatie Methode

Tegenwoordig bestaat er een groot aantal verschillende configuraties voor min kunnen we elke configuratie herleiden tot het schematische diagram in fig. 2.6. Er is één actieve elektrode en de stroom wordt gedissipeerd aarding die aangebracht wordt op de buitenzijde van de lever. Voor de verwijzen we naar subsectie 2.1.1.1.

RFA. Desalniettezoals weergegeven ter hoogte van de werking van RFA

Figuur 2.6: Schematische voorstelling van de opstelling gebruikt in huidige RFA-procedures.

De verschilpunten tussen de verschillende configuraties zijn vooral gelegen in de aard van de elektroden. Standaard eindigt de elektrode in één tip, er zijn echter ook zogeheten multitineelektrodes. Deze elektroden eindigen op verschillende tippen in een parapluvorm. Dit heeft als gevolg dat er een grotere elektrodeoppervlakte is en er een stijging plaatsvindt van het weefselvolume dat schade ondervindt [36, 35]. Daarnaast hebben veel elektroden een interne koeling via de circulatie van water of een zoutoplossing. De tip wordt met andere woorden gekoeld tot een temperatuur lager dan 25 ◦ C om het toeschroeien van het weefsel naast de tip te voorkomen [35]. Door de koeling wordt het weefsel naast de elektrodetip ook gekoeld. Dit heeft als belangrijkste gevolg dat de locatie van de maximumtemperatuur in theorie wordt opgeschoven, wat op zijn beurt resulteert in een grotere ablatieregio [36]. Om de naald op de juiste positie te brengen gebruikt men echografie. Het doel is om de tumor te vernietigen. Men opteert er vaak voor om 1 cm gezond weefsel rond de tumorsite te behandelen om zeker te zijn al het kankerweefsel vernietigd te hebben [36, 8, 35]. Bij éénmalige toediening kan RFA een ablatie veroorzaken in een gebied met een diameter van 3 cm. Bij tumoren die groter zijn dan 3 cm, zijn er meerdere toedieningen nodig. Tumoren die groter zijn dan 5 cm zijn moeilijk te behandelen met de huidige apparatuur. Een relatief grootschalige studie in Japan toont aan dat het aantal elektrode-inserties die nodig zijn bij een RFA-procedure toeneemt met de tumorgrootte [44]. De resultaten van deze studie worden samengevat in tabel 2.2.

21

Hoofdstuk 2. Radiofrequente ablatie voor levertumoren Tabel 2.2: Gemiddeld aantal elektrode-inserties in functie van de tumorgrootte [44].

Tumorgrootte [cm]

Aantal pati¨ enten (totaal = 664)

Gemiddeld aantal elektrode-inserties

0.1-2 2.1-5 >5

200 424 40

1.5 3.2 11.7

De tijdsduur van één enkele ablatie hangt opnieuw af van de gevolgde configuratie en de gewenste grootte van de ablatiezone. Typisch voert men de RFA-procedure uit gedurende 12 ` a 25 minuten [36]. Voor multitine-elektroden is een stapsgewijze expansie aangewezen: ablatie van 2 cm gedurende 5 minuten, uitbreiden naar 3 cm gedurende 5 minuten, etc. Hierdoor zouden de ablatiezones mooi in elkaar vloeien. 2.3.1.2

Complicaties

In de literatuur vindt men een mortaliteit van 0 tot 2% voor RFA bij de behandeling van levermetastasen en van 0.3 tot 0.8% bij de behandeling van HCC [29]. Majeure complicaties treden op in 6 tot 9% van de gevallen bij levermetastasen [50] en van 2.2 tot 11% bij HCC. Een studie bij 312 patiënten met hepatische tumoren die 350 RFA-procedures ondergingen, toonde 5 overlijdens aan. Deze overlijdens waren het gevolg van leverfalen, colonperforatie en driemaal een trombose van de poortader (vooral bij levercirrose). Verder waren er 37 niet-fatale ernstige complicaties, wat een incidentie oplevert van 10.6%. De complicaties bestonden uit leverabcedatie, pleuravochtuitstorting, hypoxemie, pneumothorax, subcapsulair hematoom, acute nierinsufficiëntie, hemoperitoneum en tumoruitzaaiing langs de naaldbaan [14]. Verschillende problemen ontstaan doordat de elektrode te dicht bij het diafragma of de intra-abdominale organen geplaatst werd. Dit kan vermeden worden door RFA toe te passen via open of laparoscopische weg. Er wordt in een studie bij 36% van de patiënten ook een tijdelijk postablatiesyndroom beschreven met koorts, onwel zijn, rillingen, nausea en gestegen levertransaminasen [17]. 2.3.1.3

Resultaten

In de literatuur is er een grote verscheidenheid aangaande de resultaten van RFA in de monopolaire configuratie. Algemeen kunnen we wel zeggen dat succes afhangt van enerzijds de lokalisatie van de tumor binnen de lever en anderzijds de tumorgrootte. Qua lokalisatie van de tumor speelt de nabijheid van een bloedvat een grote rol. Tumoren dicht bij grote bloedvaten worden minder succesvol behandeld door RFA. De verklaring is te vinden in een heat sink effect van de snelle bloedstroom die de warmte wegvoert van de letsels [36]. Wanneer een bloedvat als groot wordt beschouwd varieert van studie tot studie. Eenzelfde auteur gebruikt als definitie voor groot bloedvat een bloedvat met een diameter groter dan 1 cm bij een studie van RFA bij levermetastasen [31] en een diameter groter dan 3 mm bij een studie van RFA bij HCC [32]. Deze laatste studie toont een duidelijk omgekeerd verband tussen tumornecrose en de aanwezigheid van een groot bloedvat:


22

Tabel 2.3: Verband tussen tumornecrose en de aanwezigheid van een groot bloedvat [32].

Gebied

Complete tumornecrose

Perivasculair Niet-perivasculair

47% 88%

Andere probleemlocaties zijn letsels gelegen dicht bij de koepel of langs de onderrand van de lever. Deze tumoren zijn minder aangewezen voor percutane behandeling gelet op het risico voor beschadiging van het middenrif of darmperforatie. Ook de tumorgrootte speelt uiteraard een belangrijke rol. Tabel 2.4 geeft de invloed van de tumorgrootte weer in één van de vele artikels [51]: Tabel 2.4: Verband tussen tumornecrose en de tumorgrootte [51].

Tumorgrootte [cm]

Complete tumornecrose

<3 3-5 >5

80-90% 50-70% lager

We kunnen besluiten dat RFA kan toegepast worden op ambulante wijze en de meest aangewezen methode is voor patiënten met levertumoren die niet in aanmerking komen voor curatieve operatie omwille van tumorgrootte, locatie, multifocaliteit of onvoldoende leverreserve. De beste resultaten met RFA worden behaald bij patiënten met drie of minder letsels met een maximale diameter van 3.5 cm en die niet gelegen zijn dicht bij grote bloedvaten. RFA kan ook gebruikt worden bij tumorrecidief [40] en om het aantal dropouts op de levertransplantlijst te verminderen door de tumorprogressie te vertragen.

2.3.2

Nieuwe configuratie

In deze masterproef wensen we de conventionele configuratie te verlaten en onze aandacht te vestigen op de zogeheten bipolaire configuratie (fig. 2.7). In deze configuratie werken we met twee equivalente elektroden die door een stroomkring gestuurd worden. In plaats van de aarding is het nu de tweede elektrode die het elektrisch circuit sluit. Als gevolg hiervan zijn er twee ablatiezones. In tegenstelling tot de monopolaire configuratie kunnen we nu wel spelen met de stroomvormen van beide elektroden om zo een optimale overlap te bekomen tussen beide ablatiezones. In deze masterproef gaan we op zoek naar een strategie om een dergelijke overlap te bekomen. Dit zou behandeling van grotere tumoren mogelijk moeten maken.


Figuur 2.7: Schematische voorstelling van een nieuwe configuratie voor RFA.

23

Hoofdstuk 3

Numerieke methoden voor radiofrequente ablatie It is going to be necessary that everything that happens in a finite volume of space and time would have to be analyzable with a finite number of logical operations. The present theory of physics is not that way, apparently. It allows space to go down into infinitesimal distances, wavelengths to get infinitely great, terms to be summed in infinite order, and so forth; and therefore, if this proposition [that physics is computer-simulatable] is right, physical law is wrong. Richard Feynman

3.1

Inleiding

Tegenwoordig zijn heel wat commerciële softwarepakketten beschikbaar met tal van ingebouwde functies zodat het simuleren van ingewikkelde fysische fenomenen op relatief eenvoudige wijze kan plaatsgrijpen. In de literatuur omtrent numerieke simulaties van RFA wordt er frequent gebruik gemaakt van commerciële eindige elementen software of andere computationele pakketten. Wij opteren in deze masterproef om hier geen gebruik van te maken. We schrijven eigen code in C++ en de resulterende data verwerken we met eigen MATLAB-scripts. Ondanks het feit dat er veel tijd vereist is om eigen code te ontwikkelen, zijn er heel wat voordelen aan verbonden. Vooreerst is eigen software heel flexibel. Zo kunnen er onmiddellijk Magnetic Resonance Imaging (MRI)-beelden ingelezen worden zonder conversieproblemen en kunnen we de naaldposities variëren binnen eenzelfde routine om een directe vergelijking mogelijk te maken, etc. Het is verder ook mogelijk de verschillende deelproblemen te behandelen binnen eenzelfde programma. Ook zijn de koppelingen tussen de verschillende problemen direct toegankelijk. Zoals verder zal blijken, kunnen we een modulatie van de bloedperfusie met de overlevingsgraad beschouwen (twee verschillende deelproblemen). Dit is niet mogelijk met de veelgebruikte commerciële software. Ten derde hebben we ook direct zicht op de invloed van de verschillende parameters op de temperatuursdistributie binnen de tumor, wat 24

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie

25

een grote hulp is bij het opstellen van een optimale therapie. Naast een verhoogde functionaliteit is eigen code ook efficiënter. Indien we gebruik maken van geavanceerde algoritmen kunnen we de rekentijd minimaliseren en zo een verhoogde tijdsefficiëntie bekomen ten opzichte van commerciële pakketten. Dit betekent dat we in eenzelfde tijdsspanne een grotere discretisatie kunnen uitvoeren en dus nauwkeurigere resultaten verkrijgen. Ook is het zo dat we de optimale interne parameters kunnen bepalen voor onze solvers, wat de stabiliteit en numerieke nauwkeurigheid ten goede komen. Bij commerciële pakketten heb je vaak geen inzage in de gebruikte algoritmes en heb je absoluut geen invloed op interne parameters. Eigen ontwikkelde code is noodzakelijk om een nauwkeurig model op te stellen voor de nieuwe naaldconfiguratie zoals beschreven in sectie 2.3.2. De wiskundige modellering van radiofrequente ablatie in de nieuwe configuratie bestaat uit heel wat onderdelen. Zo is er een elektromagnetisch deelprobleem, een thermisch deelprobleem en een thermisch schademodel. Deze deelproblemen zijn bovendien sterk verbonden met elkaar. In dit hoofdstuk zullen we deze deelproblemen en de geassocieerde numerieke methodes belichten. Doel is om na te gaan of de nieuwe configuratie al dan niet voordelig is voor RFA.

3.2

Discretisatie van simulatieomgeving

Volgens de huidige inzichten in de natuurkunde zijn ruimte en tijd continue grootheden. Bovendien uiten de meeste andere fysische grootheden zich op macroscopische schaal als continue variabelen. Hedendaagse rekenmodules beschikken echter over een eindig geheugen en moeten methodes hanteren die in een eindig aantal logische operaties tot het resultaat komen. Daarom is de eerste stap van een numerieke simulatie de discretisatie: continue modellen, vergelijkingen en variabelen dienen omgezet te worden in hun discrete tegenhangers. Een belangrijke stap in deze discretisatie is de opdeling van de simulatieomgeving in zijn bouwstenen. In deze masterproef opteren we ervoor om het volume op te delen in kleine kubische volumes, voxels genaamd. We beschouwen deze voxels als elementaire bouwstenen, er is geen substructuur aanwezig. Elke fysische eigenschap zoals bijvoorbeeld de temperatuur in de voxel, de elektrische geleidbaarheid van de voxel, etc. wordt bijgevolg constant ondersteld over dit volume. De discretisatie is uiteraard een benadering, want in principe kunnen de fysische grootheden continu variëren binnen de voxel. Hoe kleiner het volume van de voxel, hoe accurater het eindresultaat zal zijn. Standaard werken we in deze masterproef met kubische voxels waarvoor de afmetingen ∆x = ∆y = ∆z = 0.001 m bedragen. Elke voxel kan gekarakteriseerd aan de hand van zijn relevante fysische grootheden. Aangezien cellen van hetzelfde weefsel dezelfde eigenschappen vertonen, verkiezen we om elke voxel een voxellabel te geven dat uniek is voor elk weefsel. De geometrie wordt dus opgebouwd met een lijst van cijfers, de voxellabels, die verwijzen naar de aard van het weefsel. Daarnaast wordt er een matrix opgesteld waarin voor elk weefsel de relevante fysische constanten zijn opgegeven. Een overzicht wordt gegeven in tabel 3.1. De gegevens voor gezonde en maligne levercellen zijn in overeenstemming met de waarden van tabel 2.1. Verder karakteriseren we de bloedperfusie in de lever via de fysische constanten ω, ρb en cb . Deze worden weergegeven in tabel 3.2 [8].

26


Tabel 3.1: Materiaaleigenschappen van de weefsels betrokken in de numerieke simulatie van RFA.

Type weefsel

lucht

gezond weefsel

maligne weefsel

Voxellabel σ0 S m−1 σ1 K −1 k0 W m−1 K −1 k1 W m−1 K −2 c J kg −1 K −1 ρ kg m−3

1 10−7 0 0.025 0 1012 1

2 0.36 0.015 0.53 0.00116 3500 1016

3 0.45 0.015 0.53 0.00116 3500 1016

Tabel 3.2: Fysische constanten die de bloedperfusieterm karakteriseren [8].

Grootheid ω s−1 cb J kg −1 K −1 ρb kg m−3

Numerieke waarde 0.0064 4180 1000

De naalden die de stroom injecteren worden in eerste instantie niet opgenomen in de lijst van voxellabels. We willen immers tijdens simulaties de naalden kunnen roteren om zo een optimale opstelling te verkrijgen. Daartoe moeten we de achtergrondgeometrie als basis gebruiken en voor elke positie van de naalden de achtergrond samenvoegen met de naalden. Op basis van de lijst met voxellabels kunnen we met behulp van MATLAB de simulatieomgeving op relatief eenvoudige wijze visualiseren. Een voorbeeld van een visualisatie is weergegeven in fig. 3.1. In deze figuur wordt een dwarsdoorsnede in de z-richting weerge100

y [voxelnummers]

80

60

40

20

20

40

60

80

100

x [voxelnummers]

Figuur 3.1: Dwarsdoorsnede van een simulatieomgeving.

geven. Op de figuur vinden we een doorsnede van een sferische tumor omgeven door een circulair gebied gezond weefsel. Dit geheel is omgeven door lucht. Ook de naalden zijn duidelijk merkbaar in de figuur, waarbij de tip van de naald een ander label krijgt. In de tip van de naald wordt immers de stroom ge¨ınjecteerd of geëxtraheerd, waardoor de eigenschappen


27

van de tip niet dezelfde zijn als in de rest van de naald. In principe kunnen we met herkenningssoftware zoals SPM8 patiëntspecifieke MRI-beelden inlezen. Deze software maakt het mogelijk om lever, beenderen, bloedvaten en andere weefsels te gaan differentiëren en om te zetten naar een lijst van voxels elk met het label corresponderend met het vermoedelijke weefsel [21]. In deze masterproef maken we echter vooral gebruik van een kunstmatige achtergrondgeometrie door de lever en tumoren aan de hand van wiskundige vormen te beschrijven zoals weergegeven in bovenstaande fig. 3.1. Dit verandert echter niets aan de algemeenheid van de simulatie.

3.3 3.3.1

Elektromagnetisch deelprobleem Wiskundige modellering

Het elektromagnetisch deelprobleem bestaat er in om het elektrisch veld in het doelvolume te berekenen op basis van de stroominjectie of -extractie aan de tip van beide naalden. De wiskundige vergelijking die een grote rol speelt in dit probleem is de continu¨ıteitsvergelijking (het behoud van vierstroom): ∂ρ + ∇ · J~ = 0 (3.1) ∂t Vgl. (3.1) kan eenvoudig afgeleid worden op basis van de wetten van Maxwell. Ze kan echter ook bepaald worden aan de hand van een continu¨ ummodel. Beschouwen we de netto-stroom I ([A]) in het volume V : ZZ ZZZ ~ I=− ~n · J dS = − ∇ · J~ d~v (3.2) V

∂V

Hierbij is J~ de elektrische stroomdichtheid die wordt uitgedrukt in mA2 . Behoud van lading vereist dat de netto-stroom die in het volume ge¨ınjecteerd wordt gelijk is aan de nettoladingsverandering binnen het volume ZZZ dq =I=− ∇ · J~ d~v (3.3) dt V

met q de lading ([C]) binnen het volume V of dus: ZZZ q= ρ d~v

(3.4)

V

Combinatie van vgl. (3.3) en (3.4) levert: ZZZ ∂ρ ~ + ∇ · J d~v = 0 ∂t

(3.5)

V

of er geldt in het volume V dat:

∂ρ + ∇ · J~ = 0 ∂t Om tot het elektrisch veld te komen, maken we gebruik van de wet van Ohm: ~ J~ = σ E

(3.6)

(3.7)


28

S ~ de met J~ de elektrische stroomdichtheid mA2 , σ de elektrische geleidbaarheid m en E V elektrische veldsterkte m . Met deze constitutieve wet en de definitie van de elektrische ~ = −∇φ) verkrijgen we de volgende Poisson-vergelijking: potentiaal (E ∇ · [σ (~r) ∇φ (~r, t)] =

∂ρ (~r, t) ∂t

(3.8)

Als randvoorwaarde voor deze partiële differentiaalvergelijking gebruiken we de Dirichletrandvoorwaarde: op de rand van het beschouwde volume is de potentiaal nul i.e. φ (~r, t) = 0 als ~r ∈ ∂V . Aangezien de fysische werkelijkheid zodanig is dat de potentiaal nul bedraagt op ∞, zal de Dirichlet-randvoorwaarde steeds nauwkeurigere resultaten opleveren naargelang het simulatievolume groter wordt. Het is dus aan te raden om telkens enkele lagen lucht te voorzien buiten de gesimuleerde lever om zo goed mogelijk de fysische werkelijkheid te benaderen. Bij toepassing van RFA in de nieuwe configuratie is er geen verandering van ladingsdichtheid in het beschouwde volume, tenzij op de uiteinden van de twee naalden waar de injectie of extractie van stroom plaatsgrijpt. Daarom kunnen we vgl. (3.8) herschrijven als volgt: ∇ · [σ (~r) ∇φ (~r, t)] = I0 δ (~r − ~r1 ) − I0 δ (~r − ~r2 )

(3.9)

met ~r1 en ~r2 de plaatsvectoren van de uiteinden van de twee naalden en I0 een volumetrische stroomdichtheid mA3 . Deze I0 komt overeen met de verhouding van enerzijds de stroom I doorheen de voxel die het uiteinde van de naald betreft en anderzijds het volume Vvoxel van de voxel: RRR RRR ρ d~v ρ d~v d hρivoxel d V I (3.3) 1 dq (3.4) d V RRR = = = = (3.10) I0 = Vvoxel Vvoxel dt dt Vvoxel dt dt d~v V

We onderstellen hierbij met andere woorden dat de naald een uniforme ladingsdichtheid zal induceren in de voxel waar het uiteinde zich bevindt. Het is verder ook duidelijk dat de oplossingen van vgl. (3.9) lineair zijn. Als I0 → a · I0 met a ∈ R, dan zal immers ook de oplossing eenzelfde gedrag vertonen: φ → a · φ. Eenvoudige substitutie van a · I0 en a · φ levert inderdaad dezelfde vergelijking op als de oorspronkelijke vergelijking: a · φ is dus een oplossing. Om te bewijzen dat a · φ de oplossing is moeten we gebruik maken van de uniciteitsstelling van Poisson. Stel dat zowel φ1 als φ2 oplossingen zijn van vgl. (3.9), beschouw dan φ = φ2 − φ1 . Dan voldoet φ aan: ∇ · [σ (~r) ∇φ (~r, t)] = 0

(3.11)

∇ · (φ σ∇φ) = φ∇ · (σ∇φ) + ∇φ · σ∇φ

(3.12)

Nu is

Wegens vgl. (3.11) is de eerste term van het rechterlid van de vorige vergelijking nul en dus: ∇ · (φ σ∇φ) = σ (∇φ)2 of

ZZZ V

∇ · (φ σ∇φ) d~v =

ZZZ V

σ (∇φ)2 d~v

(3.13)

(3.14)

29

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie Gebruik makend van de divergentiestelling van Gauss-Ostrogradsky verkrijgen we: ZZ ZZZ σ (∇φ)2 d~v ~n · (φ σ∇φ) dS =

(3.15)

V

∂V

Aangezien zowel φ1 als φ2 aan dezelfde Dirichlet-randvoorwaarde voldoet, zal φ = 0 op de rand van het volume. Hierdoor zal de oppervlakteintegraal over de rand van het volume in het linkerlid van vgl. (3.15) eveneens nul worden. We verkrijgen dat: ZZZ σ (∇φ)2 d~v = 0 (3.16) V

Nu is σ > 0 zodat ∇φ = 0 of ∇φ1 = ∇φ2 ↔ φ1 = φ2 + c. Aangezien φ1 = φ2 = 0 op ∂V moet φ1 = φ2 . We mogen dus inderdaad concluderen dat a · φ de unieke oplossing is als we een stroomdichtheid a · I0 aanleggen.

3.3.2

Numerieke methoden voor het elektromagnetisch deelprobleem

De Poisson-vergelijking zoals geformuleerd in vgl. (3.9) werd opgelost met behulp van de eindige differentiemethode. Deze methode gebruikt een regulier kubisch discretisatierooster als oplossingsdomein. Op die manier blijft het mogelijk om gesegmenteerde magnetische resonantie beelden te gebruiken. Op de knooppunten van dit kubisch rooster worden de afleidingsoperatoren benaderd door eindige differentieoperatoren. 2 Beschouwen we bij wijze van voorbeeld een ééndimensionale functie V . We wensen ∂∂xV2 numeriek te benaderen in het domein Ω = [0, dx ]. We gebruiken hiervoor een rooster van n equidistante punten met spatiëring hx : Ω0 = [0, hx , 2hx , · · · , dx = (n − 1) hx ] zodat xi = i·hx . Om de notatie te verlichten, schrijven we nu V (xi ) = V (i · hx ) = φi wat de functie V voorstelt ter hoogte van het roosterpunt xi . In de onderstelling dat V onbeperkt continu differentieerbaar is in Ω (V ∈ C ∞ (Ω)), kunnen we V ontwikkelen in een Taylorreeks rond het punt xi : h2x ∂ 2 V h3x ∂ 3 V h4x ∂ 4 V ∂V + + + + ··· (3.17) φi+1 = φi + hx ∂x i 2 ∂x2 i 6 ∂x3 i 24 ∂x4 i ∂V h2x ∂ 2 V h3x ∂ 3 V h4x ∂ 4 V + − + + ··· (3.18) φi−1 = φi − hx ∂x i 2 ∂x2 i 6 ∂x3 i 24 ∂x4 i

waarbij het subscript i een evaluatie voorstelt in knooppunt i met xi als ordinaat. De tweede ∂2V afgeleide ∂x2 kan dan benaderd worden door vgl. (3.17) en (3.18) op te tellen: i ∂ 2 V 1 h2x ∂ 4 V = (φi+1 − 2φi + φi−1 ) − + ··· (3.19) ∂x2 i h2x 12 ∂x4 i 2 4 Als we een truncation error toelaten van h12x ∂∂xV4 , dan kunnen we de numerieke benadering i van de tweede afgeleide schrijven als: 1 ∂ 2 V ≈ 2 (φi+1 − 2φi + φi−1 ) (3.20) 2 ∂x i hx Door vgl. (3.18) af te trekken van vgl. (3.17), verkrijgen we een centrale benadering voor de eerste afgeleide: ∂V φi+1 − φi−1 h2x ∂ 3 V = − − ··· (3.21) ∂x i 2hx 6 ∂x3 i

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie of dus:

∂V φi+1 − φi−1 ≈ ∂x i 2hx

30

(3.22)

De eerste afgeleide wordt aldus benaderd door de richtingscoëfficiënt van de rechte tussen de twee naburende punten van (xi , φi ). Voor gladde functies is deze centrale differentie de meest nauwkeurige [42]. Keren we nu terug naar vgl. (3.9), dan kunnen we deze uitdrukking uitwerken tot: ∇ · [σ (~r) ∇φ (~r, t)]

(3.9)

= =

I0 δ (~r − r~1 ) − I0 δ (~r − r~2 ) 2 ∂ φ ∂2φ ∂2φ + 2 + 2 σ (x, y, z) ∂x2 ∂y ∂z ∂σ (x, y, z) ∂φ ∂σ (x, y, z) ∂φ ∂σ (x, y, z) ∂φ + + + (3.23) ∂x ∂x ∂y ∂y ∂z ∂z

Willen we deze vergelijking uitrekenen in een knooppunt P0 van een regulier kubisch rooster, dan kunnen we de eerste en tweede afgeleiden in dit knooppunt benaderen door middel van analoge uitdrukkingen van vgl. (3.20) en (3.22). Deze eindige differentiebenaderingen worden steeds uitgedrukt in functie van de elektrische geleidbaarheden en potentialen ter hoogte van P0 en zijn zes dichtste naburen Pi (i = 1, · · · , 6). Substitutie van al deze uitdrukkingen in vgl. (3.23) levert de eindige differentievoorstelling op van de Laplaciaan in knooppunt P0 . Deze kan in compacte vorm geschreven worden als [41]: ! 6 6 X X ∇ · [σ (P0 ) ∇φ (P0 )] ≈ Ai φ (Pi ) − Ai φ (P0 ) (3.24) i=1

i=1

Voor de exacte uitdrukking van de verschillende coëfficiënten Ai verwijzen we naar [41]. We merken op dat dit oplossingsschema grote verwantschap vertoont met een elektrisch netwerk van weerstanden. In die context is vgl. (3.24) een duidelijke uiting van de stroomwet van Kirchhoff. Hierbij vervullen de coëfficiënten Ai de rol van de geleidbaarheden. Uit de bovenstaande beschouwing is het duidelijk dat de oorspronkelijke partiële differentiaalvergelijking omgevormd is tot een algebra¨ısch stelsel van vergelijkingen: [A] · [Φ] = [B]

(3.25)

Hierbij zijn Φ de onbekende potentialen, A is de systeemmatrix en in die hoedanigheid afhankelijk van de geleidbaarheden σ en de discretisatielengtes. B omvat de plaats en de groottes van de elektrische bronnen. Aangezien de benaderingen ter hoogte van een knooppunt steeds een functie zijn van de grootheden in dit knooppunt zelf en zijn zes dichtste naburen, zal de systeemmatrix A heel ijl zijn. Dit laat toe om gespecialiseerde numerieke methoden te hanteren om dit algebra¨ısch stelsel op te lossen. In deze masterproef maken we gebruik van de Preconditioned Biconjugate Gradient Stabilized Method (BiCGSTAB). Numerieke code omtrent deze Preconditioned BiCGSTAB was beschikbaar in de onderzoeksgroep EELAB, zie bijvoorbeeld [15]. Deze numerieke methode werd gebruikt vanwege zijn grote computationele efficiëntie zodat de rekentijd binnen de perken blijft.


3.4 3.4.1

31

Thermisch deelprobleem Wiskundige modellering

Bij het thermisch deelprobleem wordt er getracht om de temperatuur op elke plaats en op elk tijdstip numeriek te berekenen. Het wiskundige model dat de biologische warmteoverdracht beschrijft is ge¨ınspireerd op de Pennes bio-warmtevergelijking zoals ge¨ıntroduceerd in subsectie 2.1.2.1. Vergelijking (2.5) beschrijft echter de biologische warmteoverdracht in de lever zonder dat er RFA toegepast wordt. De invloed van een RFA-procedure zit hem in een additionele term aan de Pennes bio-warmtevergelijking. Deze term is niets anders dan 2 ~ = σ E ~ = σ |∇φ|2 W3 . We verkrijgen aldus de de elektrische vermogensdensiteit J~ · E m

volgende partiële differentiaalvergelijking: ρ (~r) c (~r)

∂T (~r, t) ∂t

= ∇ · [k (~r, T (~r, t)) ∇T (~r, t)] − ωρb cb (T (~r, t) − Tc ) 2 ~ + σ (~r, T (~r, t)) E (~r, t)

(3.26)

De term in het linkerlid beschrijft de temperatuursverandering op elke plaats en elk tijdstip. De eerste term van het rechterlid beschrijft de conductie van warmte doorheen het beschouwde volume. De tweede term is een maat voor de bloedperfusie doorheen het biologische weefsel. Deze term heeft een bloed- en leverspecifieke voorfactor die niet afhangt van plaats en tijd. We veronderstellen met andere woorden een homogene doorbloeding in het weefsel: in elke voxel is de warmteafvoerende capaciteit van de bloedstroom even groot. Dit is uiteraard niet voor elk orgaan een nauwkeurige benadering. Voor de lever kan men daarentegen wel stellen dat er een quasi homogene doorbloeding is, waardoor vgl. (3.26) geoorloofd is. De derde en laatste term van het rechterlid beschrijft het elektrisch vermogen dat aangelegd wordt tijdens een RFA-procedure. Deze term met het elektrisch veld als element is het resultaat van het elektrisch deelprobleem. Er bestaat dus een koppeling tussen het elektrisch en het thermisch probleem. Uiteraard zijn de approximaties die inherent zijn aan het model van Pennes nog steeds van toepassing (zie subsectie 2.1.2.1). Ook worden grote bloedvaten met hun heat sink effect niet correct gemodelleerd. Een oplossing voor dit probleem wordt niet teruggevonden in de literatuur. Normaliter is de enige randvoorwaarde dat de temperatuur op het begintijdstip de lichaamstemperatuur is, nl. 37 ◦ C = 310.15 K. Wij opteren er ook voor om de naalden te koelen. Dit doen we door een harde randvoorwaarde in te bouwen opdat de temperatuur in de naald 300 K zou bedragen.

3.4.2 3.4.2.1

Modified midpoint-methode Beschrijving

We opteren er voor om deze tijdsafhankelijke partiële differentiaalvergelijking op te lossen aan de hand van de modified midpoint-methode. Dit is een variant van de standaard voorwaartse Euler-methode. In concreto kan men de modified midpoint-methode onder andere gebruiken

32

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie om differentiaalvergelijkingen van het type: dT = f (t) dt

(3.27)

op te lossen. De modified midpoint-methode vertrekt van de oplossing op een vorig tijdstip en rekent de temperatuur uit op het volgende tijdstip, een tijdsperiode H na het vorige tijdstip: tvorig → tvolgend = tvorig + H

(3.28)

We verdelen dit tijdsverschil H in N equidistante delen op tijdstippen: [t0 = tvorig

t1

...

tN −1

tN = tvolgend ]

(3.29)

H . De tijdsperiode tussen de equidistante tijdstippen is dan ∆t = N De eerste stap van de modified midpoint-methode is meteen de minst nauwkeurige stap:

T (t1 ) − T (t0 ) = f (t0 ) ⇐⇒ T (t1 ) = T (t0 ) + ∆t f (t0 ) ∆t Kijken we naar de definitie van een afgeleide: d t0 + t1 t0 + t1 T (t1 ) − T (t0 ) = T =f 6= f (t0 ) lim ∆t→0 ∆t dt 2 2

(3.30)

(3.31)

dan zien we dat de modified midpoint-methode zelfs bij heel kleine ∆t een belangrijke fout maakt. Bij de volgende tijdstappen maakt ze deze niet: T (tn ) − T (tn−2 ) = f (tn−1 ) ⇐⇒ T (tn ) = T (tn−2 ) + 2∆t f (tn−1 ) 2∆t

∀n ∈ {2, . . . , N } (3.32)

In de limiet voor ∆t → 0 is deze betrekking zelfs exact. Het grote verschil met de voorwaartse Euler-methode is het feit dat de oplossing van een tijdstip niet enkel afhangt van de oplossing op de vorige tussenstap maar ook van de oplossing twee stappen terug. De modified midpointmethode is met ander woorden een tweestapsmethode. Door het feit dat de waarden op tn+1 gekoppeld zijn aan deze op tn−1 en dus niet rechtstreeks met tn , die enkel via de afgeleide in de vergelijkingen sluipt, ontstaat er een oscillerende fout. Deze fout kan gestabiliseerd worden door een Gragg-smoothing te hanteren [39]: T (tN ) =

1 [T (tN ) + T (tN −1 ) + ∆t f (tN )] 2

(3.33)

Het grote voordeel van deze Gragg-smoothing is dat de numerieke fout uitgedrukt als reeks in ∆t enkel even machten bevat. Hierdoor vermijden we oscillerend gedrag. Deze Graggsmoothing kan trouwens veralgemeend worden naar de Bulirsch-Stoer-methode [12]. Dit is echter buiten het bestek van deze masterproef. We besluiten dat voor een goede werking van de modified midpoint-methode we dus een kleine ∆t moeten hanteren om zo accuraat mogelijk de afgeleide te benaderen en een grote N om de invloed van de minst nauwkeurige stap te minimaliseren. Standaard werken we in het vervolg van deze masterproef met N = 10. Een bijkomende moeilijkheid bij toepassing van RFA is het feit dat de functie f in het rechterlid nog ruimtelijke afgeleiden bevat. In het geval van de thermische vergelijking is dit rechterlid: 2 1 ∂T ~ =f = ∇ · (k∇T ) − ωρb cb (T − Tc ) + σ E (3.34) ∂t ρc


33

De probleemterm in het rechterlid is duidelijk de ∇ · (k∇T )-term. We modelleren deze ruimtelijke afgeleiden opnieuw door gebruik te maken van eindige differentiebenaderingen op een regulier kubisch discretisatierooster, meer bepaald de centrale differentie. Voor de eenvoud kunnen we een 1D-functie y(x) beschouwen op een interval van R dat verdeeld is in equidistante roosterpunten op een afstand hx van elkaar. De eerste afgeleide van y(x) in een punt xi kan analoog benaderd worden zoals in vgl. (3.22): yi0 ≈

yi+1 − yi−1 2hx

(3.35)

De tweede afgeleide y 00 (xi ) kan opnieuw benaderd worden zoals in vgl. (3.20): yi00 ≈

yi+1 − 2yi + yi−1 h2x

(3.36)

Passen we dit toe op de conductieterm dan verkrijgen we: ∇ · (k∇T ) = k∇2 T + ∇k · ∇T T (i + 1, j, l) − 2 T (i, j, l) + T (i − 1, j, l) ≈ k (i, j, l) ∆x2 T (i, j + 1, l) − 2 T (i, j, l) + T (i, j − 1, l) + ∆y 2 T (i, j, l + 1) − 2 T (i, j, l) + T (i, j, l − 1) + ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) T (i + 1, j, l) − T (i − 1, j, l) + · 2 ∆x 2 ∆x k (i, j + 1, l) − k (i, j − 1, l) T (i, j + 1, l) − T (i, j − 1, l) + · 2 ∆y 2 ∆y k (i, j, l + 1) − k (i, j, l − 1) T (i, j, l + 1) − T (i, j, l − 1) · (3.37) + 2 ∆z 2 ∆z Hierbij staat T (i, j, l) voor de temperatuur in het roosterpunt met als cartesische co¨ ordinaten (i · dx, j · dy, l · dz). 3.4.2.2

Numerieke fouten en stabiliteit

Een inschatting van de nauwkeurigheid van numerieke methoden is primordiaal. Verder moet de methode ook resulteren in een stabiele routine waardoor er constricties moeten komen op de stapgrootten in de verschillende benaderingen. In deze subsectie gaan we in op deze zaken. De meest gebruikte definitie van een numerieke fout is via het begrip truncation error. Truncation errors van een numerieke methode zijn de fouten die veroorzaakt worden door gebruik te maken van eenvoudige approximaties voor exacte wiskundige operatoren. De enige manier om deze fouten te ontlopen is door iets exact uit te rekenen. Desondanks probeert men deze fouten te beperken door zijn toevlucht te nemen tot betere benaderingen of een groter aantal kleinere stappen. Deze truncation errors worden onderverdeeld in enerzijds lokale truncation errors en anderzijds globale truncation errors. Zoals de naam het al doet vermoeden houden lokale truncation errors enkel rekening met de fouten die in één enkele tussenstap worden ge¨ıntroduceerd. Globale truncation errors houden daarentegen rekening met de geaccumuleerde fout op een


34

bepaald tijdstip, dus met inbegrip van de fouten ge¨ıntroduceerd op voorgaande stappen. In het geval van de modified midpoint-methode is de lokale truncation error van de orde O h3 met h de stapgrootte in de methode [12]. De globale truncation error is daarentegen van de orde O h2 [12]. Deze methode is dus een grootte-orde nauwkeuriger ten opzichte van de gewone Euler-methode die een lokale truncation error van de orde O h2 heeft en orde O (h) voor de globale truncation error [19]. Naast de numerieke fouten is enige kennis omtrent de stabiliteit van een numerieke methode essentieel om een optimale stapgrootte te bepalen voor deze methode. Hoe kleiner de tussenstappen, des te groter de computationele last en dus rekentijd. Hoe groter de tussenstappen, hoe meer kans er is dat er een instabiele situatie optreedt en toevallige fouten exponentieel groeien waardoor de numerieke oplossing inaccuraat wordt. Vanuit theoretisch standpunt zijn er heel wat zienswijzen naar de stabiliteit van numerieke methoden. Een van de meest veelzijdige hulpmiddelen bij de studie van stabiliteit van eindige differentieschema’s is de Fourier methode van Von Neumann. Deze methode staat ook gekend als de Von Neumannstabiliteitsanalyse. Het grote voordeel van deze methode is zijn eenvoud: er is geen nood om de eigenwaarden van het beginwaardenprobleem te bepalen of de matrixnorm te berekenen van de systeemmatrix. Het nadeel is dat Von Neumann-stabiliteit slechts in een aantal gelimiteerde gevallen gebruikt kan worden als nodige en voldoende voorwaarde voor stabiliteit in de zin van Lax-Richtmyer (stabiliteit als en slechts als convergentie). Desondanks verkiest men de Von Neumann-stabiliteitsanalyse boven een gedetailleerde stabiliteitsanalyse voor een veel grotere klasse van problemen. De doelstelling hierbij is niet de wiskundige rigueur, maar eerder een goede gok omtrent de nodige restricties wat betreft de stapgrootte in het eindige differentieschema. Zoals vermeld in subsectie 3.4.2.1 heeft de modified midpoint-methode een vrij gecompliceerd schema dat bestaat uit verschillende voorschriften voor verschillende tussenstappen. Bovendien is het een tweestapsmethode in de zin dat niet enkel het vorige tijdstip een rol speelt in de berekening op een volgend tijdstip. Dit bemoeilijkt de berekeningen omtrent stabiliteit. In de literatuur worden bovendien geen voorbeelden teruggevonden waarin men de stabiliteit berekent van de modified midpoint-methode. Hierdoor zien we ons genoodzaakt om de Von Neumann-stabiliteit te berekenen op de klassieke Euler-methode. Aangezien deze methode eveneens voorwaarts in de tijd is en centraal in de ruimte (FTCS) zoals in de modified midpoint-methode, kunnen we vermoeden dat de Von Neumann-stabiliteit van de Euler-methode een goede schatting geeft van deze van de modified midpoint-methode. Opnieuw benadrukken we het feit dat onze enige interesse een mogelijke restrictie qua stapgrootte is. De methode van Von Neumann bestaat uit het bekijken van de evolutie van toevallige perturbaties in numerieke schema’s. Deze perturbaties kunnen een amplificatie ondergaan wat de instabiliteit veroorzaakt. Beschouwen we T n (i, j, l) wat de computationele oplossing zonder afrondfouten op tijdstip tn en in de roostercoördinaten (i, j, l) voorstelt. Deze T n (i, j, l) voldoet aan de gediscretiseerde versie van vgl. (3.26). We bekomen de volgende algebra¨ısche vergelijking gebruik makend van de FTCS Euler-methode:


T n+1 (i, j, l) − T n (i, j, l) ρ (i, j, l) c (i, j, l) = ∆t n T (i + 1, j, l) − 2 T n (i, j, l) + T n (i − 1, j, l) k (i, j, l) ∆x2 n T (i, j + 1, l) − 2 T n (i, j, l) + T n (i, j − 1, l) + ∆y 2 T n (i, j, l + 1) − 2 T n (i, j, l) + T n (i, j, l − 1) + ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) T n (i + 1, j, l) − T n (i − 1, j, l) + · 2 ∆x 2 ∆x n k (i, j + 1, l) − k (i, j − 1, l) T (i, j + 1, l) − T n (i, j − 1, l) + · 2 ∆y 2 ∆y n k (i, j, l + 1) − k (i, j, l − 1) T (i, j, l + 1) − T n (i, j, l − 1) + · 2 ∆z 2 ∆z 2 ~ − ωρb cb (T n (i, j, l) − Tc ) + σ (i, j, l) E (i, j, l)

35

(3.38)

In elke stap van een computationeel proces treden echter afrondfouten op waardoor de computationele oplossing niet T n (i, j, l) is maar eerder Ten (i, j, l). Deze waarde bevat niet alleen de gegenereerde afrondfouten van de laatste stap, maar ook de geaccumuleerde fout van eerdere stappen. De fout tussen de verwachte oplossing op basis van de gediscretiseerde algebra¨ısche vergelijkingen en de uiteindelijk verkregen computationele oplossing is dus: n (i, j, l) = T n (i, j, l) − Ten (i, j, l)

(3.39)

Aangezien Ten (i, j, l) het resultaat is van de computationele bewerking, voldoet ze eveneens aan vgl. (3.38). Deze vergelijking is lineair en dus voldoet n (i, j, l) aan: n+1 (i, j, l) − n (i, j, l) ρ (i, j, l) c (i, j, l) = ∆t n (i + 1, j, l) − 2 n (i, j, l) + n (i − 1, j, l) k (i, j, l) ∆x2 n (i, j + 1, l) − 2 n (i, j, l) + n (i, j − 1, l) + ∆y 2 n n (i, j, l + 1) − 2 (i, j, l) + n (i, j, l − 1) + ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) n (i + 1, j, l) − n (i − 1, j, l) + · 2 ∆x 2 ∆x n k (i, j + 1, l) − k (i, j − 1, l) (i, j + 1, l) − n (i, j − 1, l) + · 2 ∆y 2 ∆y n k (i, j, l + 1) − k (i, j, l − 1) (i, j, l + 1) − n (i, j, l − 1) · + 2 ∆z 2 ∆z n − ωρb cb (i, j, l)

(3.40)

De basisidee van Von Neumann is de veronderstelling dat de initiële fout ontwikkeld kan worden in een discrete Fourierreeks: 0 (i, j, l) =

ny nz nx X X X q=1 r=1 s=1

Aqrs eikq ·i∆x eikr ·j∆y eiks ·l∆z

(3.41)


36

met kq = nxπq∆x , kr = nyπr∆y en ks = nzπs ∆z . Aangezien fouten exponentieel toe- of afnemen in de tijd mogen we aannemen dat de volgende eenvoudige uitdrukking geldt: n

(i, j, l) =

=

=

ny nz nx X X X q=1 r=1 s=1 ny nz nx X X X q=1 r=1 s=1 ny nz nx X X X q=1 r=1 s=1

Aqrs eikq ·i∆x eikr ·j∆y eiks ·l∆z · eatn Aqrs eikq ·i∆x eikr ·j∆y eiks ·l∆z · ea·n∆t Aqrs eikq ·i∆x eikr ·j∆y eiks ·l∆z · ξ n

(3.42)

met ξ = ea·∆t de amplificatiefactor. De numerieke methode zal pas stabiel zijn als de fouten niet groeien of dus als |ξ| ≤ 1. Aangezien vgl. (3.40) lineair is, volstaat het om zonder verlies van algemeenheid één enkele Fouriercomponent te beschouwen van vgl. (3.42): n (i, j, l) = ξ n eikx ·i∆x eiky ·j∆y eikz ·l∆z

(3.43)

Substitutie van vgl. (3.43) in vgl. (3.40) levert: ξ−1 ρ (i, j, l) c (i, j, l) = ∆t ikx ∆x e − 2 + e−ikx ∆x eiky ∆y − 2 + e−iky ∆y eikz ∆z − 2 + e−ikz ∆z k (i, j, l) + + ∆x2 ∆y 2 ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) eikx ∆x − e−ikx ∆x · + 2 ∆x 2 ∆x ik y k (i, j + 1, l) − k (i, j − 1, l) e ∆y − e−iky ∆y + · 2 ∆y 2 ∆y ik ∆z k (i, j, l + 1) − k (i, j, l − 1) e z − e−ikz ∆z + · 2 ∆z 2 ∆z − ωρb cb (3.44) Gebruik makend van de relatie tussen de goniometrische functies en complexe exponentiëlen, verkrijgen we: ξ−1 = ρ (i, j, l) c (i, j, l) ∆t 2 (cos (kx ∆x) − 1) 2 (cos (ky ∆y) − 1) 2 (cos (kz ∆z) − 1) k (i, j, l) + + ∆x2 ∆y 2 ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) i sin (kx ∆x) · + 2 ∆x ∆x k (i, j + 1, l) − k (i, j − 1, l) i sin (ky ∆y) + · 2 ∆y ∆y k (i, j, l + 1) − k (i, j, l − 1) i sin (kz ∆z) + · 2 ∆z ∆z − ωρb cb

(3.45)


37

ξ−1 ρ (i, j, l) c (i, j, l) = ∆t 4 4 4 2 kx 2 ky 2 kz sin sin sin k (i, j, l) − ∆x − ∆y − ∆z ∆x2 2 ∆y 2 2 ∆z 2 2 − ωρb cb k (i + 1, j, l) − k (i − 1, j, l) sin (kx ∆x) +i 2 ∆x2 k (i, j + 1, l) − k (i, j − 1, l) sin (ky ∆y) + 2 ∆y 2 k (i, j, l + 1) − k (i, j, l − 1) + sin (kz ∆z) (3.46) 2 ∆z 2 Triviale uitwerking levert het volgende voorschrift voor ξ op: ∆t 4k (i, j, l) 2 kx ξ = 1− sin ∆x ρ (i, j, l) c (i, j, l) ∆x2 2 4k (i, j, l) 4k (i, j, l) 2 ky 2 kz sin sin ∆y + ∆z + ωρb cb + ∆y 2 2 ∆z 2 2 ∆t k (i + 1, j, l) − k (i − 1, j, l) +i sin (kx ∆x) ρ (i, j, l) c (i, j, l) 2 ∆x2 k (i, j + 1, l) − k (i, j − 1, l) sin (ky ∆y) + 2 ∆y 2 k (i, j, l + 1) − k (i, j, l − 1) + sin (kz ∆z) 2 ∆z 2

(3.47)

Om van een stabiele numerieke methode te kunnen spreken dient |ξ| ≤ 1 ⇐⇒ ξ · ξ ∗ ≤ 1: ∆t 4k (i, j, l) 2 kx ξ · ξ∗ = 1 − sin ∆x ρ (i, j, l) c (i, j, l) ∆x2 2 2 4k (i, j, l) 4k (i, j, l) 2 ky 2 kz + sin ∆y + sin ∆z + ωρb cb ∆y 2 2 ∆z 2 2 2 ∆t k (i + 1, j, l) − k (i − 1, j, l) + sin (kx ∆x) ρ (i, j, l) c (i, j, l) 2 ∆x2 k (i, j + 1, l) − k (i, j − 1, l) + sin (ky ∆y) 2 ∆y 2 2 k (i, j, l + 1) − k (i, j, l − 1) + sin (kz ∆z) (3.48) 2 ∆z 2 In het homogene geval (k (i, j, l) = k ∀(i, j, l) ∈ {1, . . . , nx } × {1, . . . , ny } × {1, . . . , nz }) reduceert vgl. (3.48) zich tot: ∆t 4k (i, j, l) ∗ 2 kx ξ·ξ = 1− ∆x sin ρ (i, j, l) c (i, j, l) ∆x2 2 2 4k (i, j, l) 4k (i, j, l) 2 ky 2 kz + sin ∆y + sin ∆z + ωρb cb (3.49) ∆y 2 2 ∆z 2 2 Om aan de stabiliteitsvereisten te voldoen in het homogene geval, moet: ∆t 4k (i, j, l) 4k (i, j, l) 2 kx 2 ky ∆x + ∆y sin sin ρ (i, j, l) c (i, j, l) ∆x2 2 ∆y 2 2 4k (i, j, l) 2 kz + ∆z + ωρb cb ≤ 2 sin ∆z 2 2

(3.50)

38

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie In het slechtste geval zijn sin2 voor de tijdstap ∆t wordt dan: ∆t ≤

kx 2 ∆x

= sin2

ky 2 ∆y

= sin2

kz 2 ∆z

= 1. De voorwaarde

2ρ (i, j, l) c (i, j, l) 1 1 1 4k (i, j, l) ∆x + + + ωρb cb 2 ∆y 2 ∆z 2

(3.51)

Uiteraard geldt deze voorwaarde ∀ (i, j, l) en moeten we dus de minimale tijdstap gebruiken. In het inhomogene geval is het niet zo eenvoudig om een voorwaarde af te leiden voor ∆t. 2 kxi De belangrijkste moeilijkheid is het verschijnen van zowel sin 2 ∆xi als sin (kxi ∆xi ) in vgl. (3.48). Hierdoor wordt het onmiddellijk identificeren van het slechtste geval onmogelijk gemaakt. In het slechtste geval bedraagt ξ · ξ ∗ = 1. We herschrijven vgl. (3.48): ξ · ξ ∗ = [1 − a (kx , ky , kz ) ∆t]2 + b2 (kx , ky , kz ) (∆t)2 = 1

(3.52)

Hierbij is: a (kx , ky , kz ) =

b (kx , ky , kz ) =

1 4k (i, j, l) 4k (i, j, l) 2 kx 2 ky sin ∆x + sin ∆y ρ (i, j, l) c (i, j, l) ∆x2 2 ∆y 2 2 4k (i, j, l) kz + sin2 ∆z + ωρb cb (3.53) 2 ∆z 2 k (i + 1, j, l) − k (i − 1, j, l) 1 sin (kx ∆x) ρ (i, j, l) c (i, j, l) 2 ∆x2 k (i, j + 1, l) − k (i, j − 1, l) + sin (ky ∆y) 2 ∆y 2 k (i, j, l + 1) − k (i, j, l − 1) + sin (kz ∆z) (3.54) 2 ∆z 2

Uitwerken van vgl. (3.52) levert: a2 ∆t2 − 2a∆t + 1 + ∆t2 b2 = 1 ∆t ∆t a2 + b2 − 2a = 0

(3.55)

De relevante oplossing voor de tijdstap ∆t is dus: ∆t (kx , ky , kz ) =

2a (kx , ky , kz ) 2 a (kx , ky , kz ) + b2 (kx , ky , kz )

(3.56)

Uiteraard moet deze bovengrens voor ∆t gelden voor alle kxi . We moeten met andere woorden min (∆t (kx , ky , kz )) bepalen. Dit probleem kan beschouwd worden als een driedimensionaal extremumvraagstuk. In [16] wordt aangetoond dat als een functie f : R3 −→ R een lokaal extremum bereikt in ~r, dat dan ∇f (~r) = ~0. Punten ~r waarvoor ∇f (~r) = 0 noemen we kritische punten. Merk op dat deze voorwaarde een nodige voorwaarde is, waardoor niet in elk punt die hieraan voldoet een lokaal extremum bereikt wordt. Berekenen we de kritische punten van ∆t waarvoor: ∂∆t ∂∆t ∂∆t ∇ [∆t (kx , ky , kz )] = , , = (0, 0, 0) (3.57) ∂kx ∂ky ∂kz

39


π π π dan bekomen we na een lange berekening dat ∆x , ∆y , ∆z een kritisch punt is. Bovendien toont een doorgedreven extremumonderzoek via de hoofdminorentest [16] aan dat in dit punt een lokaal minimum wordt bereikt. π π π Het eigenaardige is dat in dit punt b ∆x , ∆y , ∆z = 0, zodat: (3.56)

min (∆t (kx , ky , kz )) =

a

2

π π π ∆x , ∆y , ∆z

(3.58)

Deze voorwaarde is identiek gelijk aan de voorwaarde in het homogene geval (vgl. (3.51)). Indien we opteren voor een kleinere tijdstap dan degene bekomen in vgl. (3.51) zouden we dus zowel in het homogene als in het inhomogene geval een stabiele methode moeten verkrijgen.

3.5 3.5.1

Interpretatie van resultaten Thermische kostenfunctie

De uiteindelijke resulterende variabele in het gekoppelde elektromagnetische en thermische probleem is de temperatuur op elke plaats en op elk tijdstip. De temperatuur is echter geen éénduidige variabele bij de bestrijding van kanker. Waar we echt ge¨ınteresseerd in zijn is het antwoord op de vraag: hoeveel tumorcellen zijn vernietigd en hoeveel goede cellen hebben we hiervoor moeten opofferen? We wensen met andere woorden een wiskundig model op te stellen van de celschade. Zoals weergegeven in subsectie 2.1.3 speelt de temperatuur een belangrijke rol in de destructie van biologisch weefsel. Gedurende de beginfase van het onderzoek naar RFA werd de temperatuur dan ook als maat voor celschade gehanteerd. Eens de temperatuur in een cel een bepaalde drempelwaarde (bv. 70 ◦ C) had overschreden werd de cel als dood beschouwd. Dit is uiteraard een vrij na¨ıeve veronderstelling aangezien we met behulp van onze alledaagse intu¨ıtie dit idee naar de prullenmand kunnen verwijzen. Het is immers zo dat de blootstellingstijd een belangrijke rol speelt. Enkele uren zonnen waarbij de temperatuur ter hoogte van de huid tot 40 ◦ C oploopt, kan immers al tweedegraadsbrandwonden veroorzaken. Dezelfde wonden worden veroorzaakt door een kortstondige blootstelling aan kokend water. Op basis van een rudimentaire fysische intu¨ıtie moeten we dus rekening houden met de temperatuur, de blootstellingsduur en de biologische parameters eigen aan het beschouwde weefsel. De betrachting van de volgende subsectie is hiervoor een wiskundig model weer te geven. 3.5.1.1

Wiskundige modellering

In de literatuur worden er twee klassen van celschademodellen gebruikt. Vooreerst is er het cumulatieve equivalente minuten CEM43 -model [5]. Dit model is gebaseerd op de observatie dat per graad de temperatuur stijgt boven 43 ◦ C de benodigde tijd halveert voor een equivalent biologisch effect. De equivalente tijd voor een blootstelling aan temperaturen groter dan 43 ◦ C gedurende een tijd ti is dan ook: teq = ti · 2T −43

(3.59)

40

Hoofdstuk 3. Numerieke methoden voor radiofrequente ablatie Integratie van deze equivalente tijd levert de cumulatieve equivalente minuten (CEM43 ): CEM43 =

tZ einde 0

t · 2T (t)−43 dt

(3.60)

Wanneer deze CEM43 een bepaalde drempelwaarde overschrijdt, wordt verondersteld dat het weefsel een thermische necrose heeft ondergaan. Een tweede en meer gebruikte methode is om de thermische schade te modelleren via een Arrheniusmodel [6, 18, 8]. Voor het vervolg van deze masterproef zullen we de thermische schade bepalen aan de hand van deze methode. In de literatuur verleent men geen aandacht aan de achtergrond van het Arrheniusmodel, nochtans maakt een dergelijke beschouwing onmiddellijk duidelijk welke inherente veronderstellingen men toepast bij het opstellen van dit wiskundig model. Hiertoe maken we eerst een uitstap naar de wereld van de scheikunde. Beschouw de volgende scheikundige reactie: aA + bB −→ C

(3.61)

dan definieert men de reactiesnelheid r als: 1 d [B] d [C] 1 d [A] =− = (3.62) a dt b dt dt met [X] de concentratie van stof X mol . Vanuit het oogpunt van de botsingtheorie is l deze reactiesnelheid r proportioneel met de aanwezige concentraties van de reactanten: r=−

r = k (T ) [A]a [B]b

(3.63)

De proportionaliteitsfactor k (T ) wordt in scheikundige middens ook wel eens de reactiesnelheidsconstante genoemd. Deze benaming als zou k (T ) een constante zijn is een beetje misleidend: deze proportionaliteitsfactor is immers nog temperatuursafhankelijk. De vergelijking van Arrhenius beschrijft deze afhankelijkheid voor een chemische reactie zonder intermediaire reacties, de zogenaamde éénstapsreactie: ∆Ea

k (T ) = Ae− R·T

(3.64)

J en R de Hierbij is A de frequentiefactor s−1 , ∆Ea de molaire activeringsenergie mol J molaire gasconstante R = 8.314472 K·mol . De idee van Arrhenius is dat de fractie moleculen die genoeg energie hebben om de reactie te ondergaan, zal afhangen van de verhouding van de gemiddelde molaire thermische energie (R · T ) en de molaire activeringsenergie ∆Ea . De reden waarom er een exponentieel verband is, wordt verklaard door de aanname dat er thermisch evenwicht heerst en de moleculen dus een Maxwell-Boltzmannverdeling volgen. De frequentiefactor A wordt in de botsingstheorie fysisch verklaard als zijnde de frequentie of dus de kans waarmee moleculen met elkaar botsen. Deze factor heeft een verwaarloosbare temperatuursafhankelijkheid en mag voor elke specifieke reactie constant verondersteld worden. Keren we nu terug van dit intermezzo dan kunnen we het proces dat een levende cel ondergaat


41

∆E

a beschrijven als de eerste orde-reactie L −−−→ D met L als levende cel en D als dode cel. Nu geldt:

r

(3.62)

=

(3.63)

=

−

d [L] dt (3.64)

∆Ea

k (T ) [L] = Ae− R·T [L]

(3.65)

Noteren we [L] nu als c (~r, t) als de concentratie van levende moleculen op tijdstip t en op plaats ~r dan verkrijgen we door het herschrijven van vgl. (3.65): ∂c (~r, t) − ∆Ea = −Ae R·T (~r,t) · c (~r, t) ∂t

(3.66)

Deze differentiaalvergelijking kunnen we herschrijven als: ∂ ln [c (~r, t)] − ∆Ea = −Ae R·T (~r,t) ∂t

(3.67)

waardoor de oplossing c (~r, t) voldoet aan:

Zt c (~r, t) − ∆Ea ln = − Ae R·T (~r,t) dt = −Ω (~r, t) c (~r, 0)

(3.68)

0

waarbij we de Arrhenius-schadeindex Ω (~r, t) definiëren als: Ω (~r, t) =

Zt

∆Ea − R·T (~ r ,t)

Ae

dt

(3.69)

0

c(~ r,t) Beschouwen we nu α (~r, t) = c(~ r,0) . Dit is de verhouding van de aanwezige concentratie levende cellen op tijdstip t ten opzichte van de initiële concentratie levende cellen op het begintijdstip. Deze verhouding drukt met andere woorden de probabiliteit uit dat een molecule op tijdstip t nog leeft. We kunnen α dus ook de overlevingsgraad noemen:

α (~r, t) = exp [−Ω (~r, t)]

(3.70)

Op t = 0 bedraagt de overlevingsgraad uiteraard 1, gedurende de thermische belasting zal deze overlevingsgraad monotoon afnemen en naderen naar 0. Aangezien een exponentiële functie geen nulpunten heeft, wordt er vaak een drempelwaarde ge¨ınstalleerd voor de overlevingsgraad. Eens α onder deze waarde komt, stellen we α = 0 en beschouwen we de cel als vernietigd. In de literatuur gebruikt men vaak de drempelwaarde αdrempel = 0.01 wat overeenkomt met een kans van 99 % dat de cellen vernietigd worden [18]. Vaak drukt men deze kans uit via de Arrhenius-schadeindex: Ω ≥ 4.6 (exp (−4.6) ≈ 0.01). Een andere frequente drempelwaarde is Ω = 1 wat dan overeenkomt met de kans van 63 % dat de cellen op de beschouwde plaats dood zijn [6, 5, 8]. Deze drempelwaarde wordt vaak in verband gebracht met het tijdstip waarop thermische necrose optreedt [6]. Desondanks verkiezen we om de drempelwaarde αdrempel = 0.01 te hanteren gezien de interpretatie en het veelvuldige gebruik in de literatuur. Er rest natuurlijk nog het probleem om de parameters A en ∆Ea te bepalen voor cellen in de lever. Net zoals de materiaaleigenschappen in subsectie 2.2.4 is het een moeilijke opdracht


42

om deze parameters accuraat in te schatten. We durven te stellen dat dit naast de wiskundige modellering van biologische systemen de grootste uitdaging is van de biomedische wetenschappen. Er zijn enorm veel invloeden op cellulair vlak en het is niet evident om experimenten uit te voeren op levende biologische systemen. In de literatuur verschijnen er dan ook tal van uiteenlopende parametersets [6, 18, 8]. Een overzicht van de twee populaire sets wordt gegeven in tabel 3.3. Tabel 3.3: Overzicht van de parameters gebruikt in het Arrheniusmodel voor thermische schade van biologische weefsels.

Bron A s−1 J ∆Ea mol

[18]

[8]

1080

3.1 · 1098 6.28 · 105

2.984 · 5.064 · 105

Wegens het veelvoudige gebruik van de set en de autoriteit van de originele auteurs, gebruiken we de parameters zoals beschreven in [18]. Enige gegevens omtrent de experimentele validatie van deze of andere parameters worden niet teruggevonden in de literatuur. 3.5.1.2

Numerieke methode

In deze subsectie beschrijven we hoe vgl. (3.70) numeriek wordt berekend op basis van de temperatuurdistributie. Bekijken we de Arrhenius-schadeindex op tn : Ω (~r, tn )

(3.69)

=

Ztn

Ae

∆Ea − R·T (~ r ,t)

dt

t0

=

ti+1 Z − ∆Ea Ae R·T (~r,t) dt

n−1 X

i=0 t i

=

n−1 X

Ωi (~r)

(3.71)

i=0

Voor de overlevingsgraad betekent dit dat: α (~r, tn )

(3.70)

=

(3.71)

=

exp [−Ω (~r, tn )] " n−1 # X exp − Ωi (~r) i=0

=

n−1 Y

exp [−Ωi (~r)]

(3.72)

i=0

Door deze eigenschap leent de overlevingsgraad zich tot een recursieve berekening: α (~r, tn )

(3.72)

=

=

α (~r, tn−1 ) · exp [−Ωn−1 (~r)]   Ztn   − ∆Ea α (~r, tn−1 ) · exp − Ae R·T (~r,t) dt tn−1

(3.73)

43


We kunnen de overlevingsgraad van de huidige stap dus berekenen op basis van de overlevingsgraad van de vorige stap, wat uiteraard computationeel heel efficiënt is. De bepalende stap is Rtn − ∆Ea dus de numerieke berekening van Ωn−1 (~r), i.e. de numerieke integratie van Ae R·T (~r,t) dt. tn−1

Voor de numerieke integratie van functies bestaan er heel wat verschillende methoden. Wij opteren voor de eenvoudige trapeziumregel omdat we toch maar beschikken over de temperatuur op ti en ti+1 . Er zijn geen tussenliggende temperaturen beschikbaar en dus is een Gaussische kwadratuur bijvoorbeeld niet zinvol. Alternatief zouden we kunnen afstappen van de recursiviteit en een Gaussische kwadratuur kunnen toepassen op de volledige integraal van t0 tot tn . Aangezien we op elk tijdstip de overlevingsgraad wensen te kennen, zou dit resulteren in een verminderde computationele efficiëntie. De algemene werking van de trapeziumregel bestaat er in de functie f (x) op het gesloten interval [a, b] te benaderen door de lineaire functie met eenzelfde eindpunten: f (x) ≈

b−x x−a f (a) + f (b) b−a b−a

(3.74)

De integratie van f (x) over [a, b] is dan: Zb a

f (x) dx

(3.74)

≈

= =

Z

a

b

x−a b−x f (a) + f (b) dx b−a b−a

b−a 1 (bf (a) − af (b)) + b−a b−a b−a (f (a) + f (b)) 2

b2 a2 − 2 2

(f (b) − f (a)) (3.75)

wat uiteraard de oppervlakte is van de benaderende trapezium, nomen est omen. Passen we dit toe op Ωn−1 (~r) dan verkrijgen we: Ωn−1 (~r) =

Ztn

tn−1

≈

A exp −

∆Ea dt R · T (~r, t)

∆t ∆Ea ∆Ea A exp − + exp − 2 R · T (~r, tn−1 ) R · T (~r, tn )

(3.76)

met ∆t = tn − tn−1 . Substitutie van vgl. (3.76) in vgl. (3.73) levert: α (~r, tn ) ≈ α (~r, tn−1 ) ∆Ea ∆Ea ∆t + exp − (3.77) · exp − A exp − 2 R · T (~r, tn−1 ) R · T (~r, tn )

3.5.2

Statistische verwerking

Eens we de overlevingsgraad kunnen bepalen, dan zijn we in staat om het resultaat van de RFA-procedure statistisch te verwerken. Uiteraard willen we bepalen hoeveel tumorcellen er procentueel verdwenen zijn. Het is echter ook interessant om eens te bekijken hoeveel gezonde cellen er vernietigd zijn. Om dergelijke zaken statistisch te verwoorden, kunnen we analoge begrippen opstellen zoals voor het uitdrukken van de diagnostische waarde van een test [4]. We hanteren de volgende begrippen: waar positief (TP), vals positief (FP), waar


44

negatief (TN) en vals negatief (FN). De definitie van deze begrippen in de context van RFA wordt weergegeven in tabel 3.4. Tabel 3.4: Definitie van de statistische parameters die de uitkomst van een RFA-procedure karakteriseren.

Tumorcel

Gezonde cel

TP FN

FP TN

Verbrand Niet verbrand

Verder definiëren we de waar positieve fractie (TPF) als de verhouding van de verbrande tumorcellen ten opzichte van alle oorspronkelijke tumorcellen: TPF =

TP TP + FN

(3.78)

In statistisch jargon staat de TPF gekend als de sensitiviteit [4]. Analoog is de vals positieve fractie (FPF) de verhouding van de verbrande gezonde cellen en de totale hoeveelheid gezonde cellen: FP FPF = (3.79) TN + FP Statistisch gezien is de specificiteit interessanter [4]: Specificiteit = 1 − F P F =

TN TN + FP

(3.80)

Dit is de verhouding van de gezonde cellen die niet aangetast zijn door RFA ten opzichte van alle oorspronkelijke gezonde cellen. Om de doelstelling van de RFA-procedure te beoordelen kunnen we gebruik maken van de statistische nauwkeurigheid: Nauwkeurigheid =

TP + TN TP + FP + TN + FN

(3.81)

Dit drukt de verhouding uit van de cellen waarbij de doelstelling bereikt is ten opzichte van alle cellen. Een alternatieve statistische parameter is de gebalanceerde nauwkeurigheid: Gebalanceerde nauwkeurigheid

= (3.80)

=

=

sensitiviteit + specificiteit 2 1 (T P F + 1 − F P F ) 2 1 TP TN + 2 TP + FN TN + FP

(3.82)

Deze parameter geeft het rekenkundig gemiddelde weer van enerzijds de fractie tumorcellen die verbrand zijn en anderzijds de fractie gezonde cellen die niet verbrand zijn. Aangezien het de doelstelling is van een RFA-procedure om zo veel mogelijk tumorcellen te doden en tegelijkertijd zo veel mogelijk gezonde cellen te sparen, lijkt de gebalanceerde nauwkeurigheid ons de uitverkoren parameter. Om te kunnen beoordelen of de ene configuratie voordeliger is ten opzichte van een andere zullen we deze gebalanceerde nauwkeurigheid hanteren.


3.6

45

Modulatie van de Pennes bio-warmtevergelijking

Zoals aangegeven in subsectie 2.1.3, neemt men waar dat de microvasculatuur lijdt onder de thermische schade. De perfusie door de microvasculatuur neemt af en indien er voldoende thermische schade aangericht is, ontstaat er een vasculaire shutdown (geen doorstroming). Het lijkt ons dan ook niet correct om gedurende RFA te werken met de bloedperfusieterm ωρb cb (T (~r, t) − Tc ) zoals beschreven in vgl. (3.26). Aangezien ω, ρb en cb plaatsen tijdsonafhankelijk zijn, onderstelt men impliciet dat de doorbloeding homogeen is in elke voxel. Nu stelt de overlevingsgraad α de verhouding voor van de aanwezige concentratie van levende cellen ten opzichte van de initiële concentratie. Indien we de zware veronderstelling maken dat de microvasculatuur eenzelfde schademodel volgt als de individuele cellen, dan kunnen we stellen dat α ook de verhouding voorstelt van werkzame bloedvaten op microschaal ten opzichte van de initiële hoeveelheid werkzame bloedvaten op microschaal. De veronderstelling dat de bloedperfusie in elke voxel homogeen is, impliceert dat er een homogene distributie is van bloedvaten op microschaal. Beschouw nu n bloedvaten. Onderstel nu dat deze een gezamenlijk warmteafvoerend effect hebben: Qaf voer . Na een tijd t is er een overlevingsgraad α (t), wat impliceert dat er slechts n · α (t) bloedvaten meer aanwezig zijn. Het warmteafvoerend effect zal daarom gereduceerd worden tot Qaf voer · α (t). Deze primitieve redenering zet ons aan tot een modulatie van de bloedperfusieterm in het model volgens Pennes met de overlevingsgraad: ωρb cb (T (~r, t) − Tc ) −→ α (~r, t) · ωρb cb (T (~r, t) − Tc )

(3.83)

Het effect van de bloedperfusieterm zal gezien de dalende overlevingsgraad steeds afnemen naarmate de thermische schade toeneemt. Aangezien de bloedvaten steeds minder warmte kunnen afvoeren, zullen we in ons model hogere temperaturen verkrijgen ten opzichte van het originele model (bij eenzelfde stroominjectie). Uiteraard zijn we ons bewust van het feit dat een simpele lineaire modulatie een proces met een dergelijke complexiteit onmogelijk accuraat kan beschrijven. We geloven echter dat dit een betere poging is dan de klassieke bloedperfusieterm die de warmteafvoerende capaciteit overschat. Door deze overschatting komt men kunstmatig lagere temperaturen uit als resultaat van de computationele berekeningen. Het is steeds gevaarlijk om de belangrijkste variabele in het fysisch proces te onderschatten. Zo kan er meer schade aangericht worden dan voorspeld door het computermodel. We vermoeden bovendien dat de overlevingsgraad van de microvasculatuur minder snel zal dalen dan deze van de omliggende cellen, gezien de beschermende structuur van de bloedvatwanden. De warmteafvoerende capaciteit zal vermoedelijk hoger zijn dan in ons model, wat op zijn beurt impliceert dat onze berekende temperaturen wellicht te hoog zijn. We overschatten met andere woorden de thermische schade. Om een accurate uitspraak te kunnen doen over deze schade hebben we nood aan een grote hoeveelheid nauwkeurige experimentele data. Zoals eerder aangegeven is dit uiteraard niet triviaal gezien de aard van het studieobject, nl. een biologisch systeem. Ondanks de logica achter dit model, leert nader onderzoek in de literatuur ons dat slechts één auteur eenzelfde modulatie voorstelt [6]. Sedert het verschijnen in 2004 wordt dit artikel niet opgevolgd in de recente literatuur. De beweegredenen zijn hiervoor niet duidelijk. Wellicht


46

willen veel wetenschappers een min of meer aanvaard model niet inruilen voor een aangepast model zonder een duidelijke bevestiging via het experiment. Bovendien is het zo dat vele artikels het computationeel model in commerciële software ontwikkelen. De verschillende deelproblemen worden hierbij meestal in verschillende programma’s gebruikt. Dit maakt het invoeren van een additionele koppeling niet evident. Het loont dus de moeite om eigen software te ontwikkelen zoals uitvoerig beargumenteerd in de inleiding van dit hoofdstuk.

3.7

Overzicht van de globale werking van het numeriek model voor radiofrequente ablatie

De dominerende wiskundige vergelijkingen van elk deelprobleem dienen samen opgelost te worden. Voor de volledigheid geven we de verschillende vergelijkingen opnieuw weer:    ∇ · [σ (~r, T (~r, t)) ∇φ (~r, t)] = I0 δ (~r − r~1 ) − I0 δ (~r − r~2 )     ∂T (~r, t)   = ∇ · [k (~r, T (~r, t)) ∇T (~r, t)] − α (~r, t) ωρb cb (T (~r, t) − Tc ) ρ (~r) c (~r)   ∂t   2 ~ (3.84) + σ (~ r , T (~ r , t)) E (~ r , t)        Zt   − ∆Ea 0  R·T (~ r ,t )   dt0  α (~ r , t) = exp − Ae    0

Zoals merkbaar in de verschillende argumenten bestaat er een sterke koppeling tussen de verschillende deelproblemen. Hierdoor ontstaat een complex geheel waarin men door de bomen het bos niet meer dreigt te zien. Daartoe geven we in fig. 3.2 een schematisch overzicht van de methode in globo.

47


Elektromagnetisch deelprobleem

I0

       T     I 0  r  r1    r  r2 

σ(T)

φ

E

Thermisch deelprobleem T

c

T    k T  T   b cb T  Tc  t 2   T  E

k(T)

Thermische kostenfunctie  t  E     exp   A exp   a  dt '  RT    0

α

Figuur 3.2: Overzicht van de globale werking van het numeriek model voor radiofrequente ablatie.

Hoofdstuk 4

Numerieke analyse van numerieke methoden Primum non nocere.

4.1

Inleiding

Aangezien RFA toegepast wordt op de lever, een belangrijk onderdeel van het menselijk lichaam, moeten we de nodige voorzichtigheid in acht nemen. We willen in ieder geval voorkomen dat de procedure het ziektepatroon erger maakt of dat er levensbedreigende complicaties zouden optreden. Numerieke methoden vormen benaderingen voor de exacte wiskundige modellen. Het is dan ook belangrijk om het resultaat van deze methoden te vergelijken met analytische uitdrukkingen. Gezien de complexiteit van de vergelijkingen is dit niet altijd eenvoudig. Desalniettemin zullen we in dit hoofdstuk voor elk afzonderlijk deelprobleem een analytische limietsituatie proberen te bedenken. Telkens zal er een kritische analyse volgen teneinde de correctheid van de implementatie te bevestigen. Als elk afzonderlijk deelprobleem correct werkt, kunnen we vermoeden dat het computermodel ons globaal wiskundig model correct uitvoert.

4.2 4.2.1

Elektromagnetisch deelprobleem Analytische limietsituatie

Het is allesbehalve eenvoudig om een elektromagnetisch probleem te bedenken waarbij er een injectie van stroom plaatsvindt en dat nog analytisch berekenbaar is. Een mogelijke limietsituatie ontstaat als we een homogene sfeer met elektrische geleidbaarheid σ in lucht plaatsen waarbij de elektroden op de rand van deze homogene sfeer stroom injecteren. 4.2.1.1

Delta-bron

We beschouwen vooreerst de onrealistische situatie waarbij we stroom injecteren of extraheren in de infinitesimaal punten die overeenstemmen met de noord- en zuidpool van de sfeer. 48

49

Hoofdstuk 4. Numerieke analyse van numerieke methoden In de sfeer, reduceert vgl. (3.9) zich tot: ∇ · (σ ∇φ (~r)) = 0

(4.1)

wat gezien de constante elektrische geleidbaarheid σ resulteert in de Laplace-vergelijking: ∇2 φ (~r) = 0

(4.2)

Door de sferische symmetrie is het aangewezen om de cartesische coördinaten ~r in te ruilen voor sferische co¨ ordinaten (r, φ, θ). Om mogelijke verwarring met de azimutale hoek φ te voorkomen noteren we de potentiaal in het vervolg van dit hoofdstuk als V . Vooreerst dienen we de Laplaciaan uit te drukken in bolcoördinaten [30]: 1 ∂ ∂ 1 ∂ ∂2 1 2 ∂ 2 r + 2 2 sin θ (4.3) ∇ ≡ 2 + r ∂r ∂r ∂θ r sin θ ∂φ2 r2 sin θ ∂θ We stellen voor om vgl. (4.2) op te lossen via een scheiding der veranderlijken: V (r, φ, θ) = R (r) Φ (φ) Θ (θ)

(4.4)

Substitutie van uitdrukkingen (4.3) en (4.4) in vgl. (4.2) geeft na deling door vgl. (4.4): 1 d 1 1 d2 Φ (φ) 1 1 d dΘ (θ) 2 dR (r) r + + sin θ =0 (4.5) R (r) dr dr sin θ Θ (θ) dθ dθ sin2 θ Φ (φ) dφ2 Aangezien de randvoorwaarde de noord- en zuidpool van de sfeer betreft, is deze randvoorwaarde onafhankelijk van de azimutale hoek φ. De oplossing van vgl. (4.5) zal dan ook niet afhangen van deze hoek: dΦ(φ) = 0 [30]. Hierdoor bekomen we de volgende partiële dφ differentiaalvergelijking: 1 1 d dΘ (θ) 1 d 2 dR (r) r =− sin θ (4.6) R (r) dr dr sin θ Θ (θ) dθ dθ Aangezien het linkerlid van deze vergelijking enkel afhangt van r en het rechterlid enkel afhankelijk is van θ moeten beide leden gelijk zijn aan een constante k. We bekomen met andere woorden twee gewone differentiaalvergelijkingen: d2 R (r) dR (r) 1 d 2 dR (r) r = k ⇐⇒ r2 + 2r − kR (r) = 0 (4.7) 2 R (r) dr dr dr dr 1 d dΘ (θ) sin θ + kΘ (θ) = 0 (4.8) sin θ dθ dθ Stellen we nu dat k = n (n + 1) met n ∈ N dan kunnen we vgl. (4.7) identificeren met een Euler-Cauchy-vergelijking. Deze differentiaalvergelijking heeft als algemene oplossing [30]: R (r) = an rn + bn r−(n+1)

met an en bn constanten

(4.9)

Vervolgens lossen we vgl. (4.8) op. Daartoe stellen we cos θ = u, dan is sin2 θ = 1 − u2 en via de kettingregel weten we dat: d d du d = = − sin θ dθ du dθ du

(4.10)

Hoofdstuk 4. Numerieke analyse van numerieke methoden Substitutie van deze betrekkingen in vgl. (4.8) levert: dΘ (u) d + n (n + 1) Θ (u) = 0 1 − u2 du du

50

(4.11)

Deze differentiaalvergelijking is niets anders dan de gekende Legendre-vergelijking [30]. Indien n ∈ N is de oplossing een Legendre-polynoom van graad n: Θ (u) = Pn (u) ⇐⇒ Θ (θ) = Pn (cos θ)

(4.12)

De algemene oplossing voor ons probleem is dus: V (r, θ) =

∞ X

n=0

an rn + bn r−(n+1) Pn (cos θ)

(4.13)

Aangezien we ge¨ınteresseerd zijn in de potentiaal binnen de sfeer en r = 0 tot dit gebied behoort, moeten we bn = 0 stellen ∀n ∈ N om een singulariteit te vermijden. De oplossing van de Laplace-vergelijking neemt voor ons probleem dus de volgende vorm aan: V (r, θ) =

∞ X

an rn Pn (cos θ)

(4.14)

n=0

Vervolgens bekijken we de randvoorwaarde. We wensen een randvoorwaarde te hanteren die stroom injecteert ter hoogte van de rand van de sfeer. Een manier om een accurate randvoorwaarde te verkrijgen is door de stroomdichtheid op te leggen ter hoogte van de volledige rand van de sfeer (r = R met R de straal van de sfeer). We gebruiken hiertoe de wet van Ohm: ~ (R, θ) = −σ∇V (R, θ) J~ (R, θ) = σ E (4.15) Scalaire vermenigvuldiging met de normaalvector ~1n = ~1r levert: J (R, θ) = −σ

∂V (R, θ) ∂r

(4.16)

Definiëren we nu J (cos θ) als de stroomdichtheid die ge¨ınjecteerd wordt (richting −~1r ) ter hoogte van de rand van de sfeer, dan komen we tot: J (cos θ) = σ

∂V (R, θ) ∂r

(4.17)

Indien we werken met een delta-bron ter hoogte van noord- en zuidpool bedraagt J (cos θ): J (cos θ) = J0 [δ (cos θ − 1) − δ (cos θ + 1)]

(4.18)

Hierbij is J0 de stroomdichtheid die ge¨ınjecteerd wordt door de elektroden. Aangezien de Laplace-vergelijking een speciaal geval is van de Poisson-vergelijking, zijn de oplossingen voor de potentiaal ook lineair met de stroom (zoals bewezen in subsectie 3.3.1). Zonder aan algemeenheid in te boeten, mogen we dus J0 = 1 stellen. We zullen in dit hoofdstuk dan ook altijd genormeerde potentialen bekijken teneinde een eerlijke vergelijking te kunnen uitvoeren. Invullen van deze betrekkingen in randvoorwaarde (4.17) levert: J (cos θ)

(4.18)

=

(4.17)

=

δ (cos θ − 1) − δ (cos θ + 1) ∞ ∂V (4.14) X (R, θ) = σ n an Rn−1 Pn (cos θ) σ ∂r n=0

(4.19)

51

Hoofdstuk 4. Numerieke analyse van numerieke methoden

Vermits J (cos θ) gelijk moet zijn aan een reeks van Legendre-polynomen, wensen we J (cos θ) zelf te ontbinden in een dergelijke reeks. We willen met andere woorden de Fourier-Legendrereeks opstellen van J (cos θ) [30]. Vermits de Legendre-polynomen een complete orthogonale set vormen, kunnen we analoog als een Fourier-reeks deze Fourier-Legendre-reeks bekomen. Stellen we: ∞ X J (cos θ) = cn Pn (cos θ) (4.20) n=0

dan kunnen we de coëfficiënten cn berekenen door linker- en rechterlid te vermenigvuldigen met Pm (cos θ) en te integreren van −1 tot 1. We verkrijgen: Z1

J (cos θ) Pm (cos θ) d cos θ =

Z1 X ∞

cn Pn (cos θ) Pm (cos θ) d cos θ

(4.21)

−1 n=0

−1

De Legendre-polynomen voldoen aan de volgende orthogonaliteitsrelatie [47]: Z1

Pn (cos θ) Pm (cos θ) d cos θ =

2 δn,m 2m + 1

(4.22)

−1

Hierdoor wordt het rechterlid van vgl. (4.21) gereduceerd tot in vgl. (4.21) levert: cm

=

2m + 1 2

Z1

2 2m+1 cm .

Substitutie hiervan

J (cos θ) Pm (cos θ) d cos θ

−1 (4.18)

=

2m + 1 (Pm (1) − Pm (−1)) 2

(4.23)

Nu geldt dat [47]: Pm (−x) = (−1)m Pm (x) Pm (1) = 1

∀x ∈ [−1, 1]

∀m ∈ N

(4.24) (4.25)

zodat we vgl. (4.23) kunnen herschrijven als: cm =

2m + 1 (1 − (−1)m ) 2

(4.26)

Vervangen we deze uitdrukking in randvoorwaarde (4.19), dan bekomen we: σ

∞ X

n an R

n−1

(4.19)

(4.26)

Pn (cos θ) = J (cos θ) =

n=0

∞ X 2n + 1

n=0

2

(1 − (−1)n ) Pn (cos θ)

(4.27)

Aangezien de Legendre-polynomen orthogonaal zijn, moet deze betrekking gelden voor elke individuele term van de reeks: σ n an Rn−1 = an =

2n + 1 (1 − (−1)n ) 2 2n + 1 1 (1 − (−1)n ) n σRn−1 2

∀n ∈ N0

(4.28)

52


Ten slotte substitueren we de uitdrukking voor deze coëfficiënten in de algemene oplossing (4.14): ∞ X 2n + 1 1 (1 − (−1)n ) n V (r, θ) = r Pn (cos θ) (4.29) n σRn−1 2 n=1

Deze oplossing kunnen we verder vereenvoudigen vermits (1 − (−1)n ) = 0 voor even n. Beschouwen we enkel de oneven termen van de machtreeks, dan bekomen we de uiteindelijke oplossing: ∞ X 4m + 3 1 V (r, θ) = r2m+1 P2m+1 (cos θ) (4.30) 2m + 1 σR2m m=0

4.2.1.2

Bolkapelektrode

De situatie met de delta-bron is onrealistisch omwille van het gegeven dat we werken met een gediscretiseerde simulatieomgeving. Een voorbeeld van een dergelijke discretisatie is gegeven in fig. 4.1.

α

R’

Figuur 4.1: Schematische voorstelling van de analytische situatie aangaande de bolkapelektrode.

Stroominjectie of extractie dient te gebeuren in minimaal één voxel (aangegeven in het groen in fig. 4.1). De verwachting is dan ook dat de delta-bron niet volledig in overeenstemming is met de numerieke simulatie. In deze subsectie zullen we trachten een analytische limietsituatie te bedenken met enige ruimtelijke uitgestrektheid. Beschouw weerom fig. 4.1. We injecteren nu stroom langsheen de booglengte van de cirkelsector met hoek α zoals aangeduid in fig. 4.1. Door rekening te houden met de azimutale symmetrie komt deze situatie overeen met de injectie van stroom doorheen een bolkap, vandaar de bolkapelektrode. Analytisch blijft naast de algemene oplossing van de Laplace-vergelijking ook de randvoor-

53

Hoofdstuk 4. Numerieke analyse van numerieke methoden waarde gelden:

 ∞ X    an rn Pn (cos θ)  V (r, θ) =

n=0 (4.31)   ∂V   J (cos θ) = σ (R, θ) ∂r Hierbij worden dezelfde definities gehanteerd als in de subsectie omtrent de delta-bron. Het enige verschilpunt met de vorige analyse betreft de vorm van J (cos θ). Er wordt nu niet enkel meer een stroomdichtheid ge¨ınjecteerd bij cos θ = ± − 1. Voor een voorstelling van deze stroomdichtheid verwijzen we naar fig. 4.2.

J/J0

1

-1 -cos α cos α 1

cos θ

-1 Figuur 4.2: Grafische weergave van het stroomdichtheidsprofiel ge¨ınjecteerd door een bolkapelektrode.

Er wordt nu stroom ge¨ınjecteerd vanaf θ = 0 tot θ = α. De grootte van deze hoek α kan 0 eenvoudig afgeleid worden op basis van fig. 4.1. We zien dat tan α = R00.5 +0.5 met R de gediscretiseerde straal van de bol zoals aangegeven in fig. 4.1. Dan is cos α: cos α = √

R0 + 0.5 =q 1 + tan2 α (R0 + 0.5)2 + 0.52 1

(4.32)

Zoals in de vorige afleiding ontwikkelen we dit stroomdichtheidsprofiel in zijn FourierLegendre-reeks: ∞ X J (cos θ) = cm Pm (cos θ) (4.33) m=0

We volgen opnieuw dezelfde strategie: Z1

J (cos θ) Pn (cos θ) d cos θ

=

Z1 X ∞

cm Pm (cos θ) Pn (cos θ) d cos θ

(4.34)

−1 m=0

−1

(4.22)

=

cn

2 2n + 1

(4.35)

54

Hoofdstuk 4. Numerieke analyse van numerieke methoden Substitutie van het profiel uit fig. 4.2 in vgl. (4.35) levert: Z1

cos α

Pn (cos θ) d cos θ −

−Zcos α

Pn (cos θ) d cos θ = cn

2 2n + 1

(4.36)

−1

Vervangen we θ in de tweede integraal door θ0 = π − θ, dan verkrijgen we: Z1

cos Z α

Pn (cos θ) d cos θ +

cos α

Pn − cos θ0 d cos θ0 = cn

1

2 2n + 1

(4.37)

Gezien vgl. (4.24) geldt er dat Pn (− cos θ0 ) = (−1)n Pn (cos θ0 ). Aldus bekomen we dat: Z1

Pn (cos θ) d cos θ +

cos α

cos Z α 1

(−1)n Pn cos θ0 d cos θ0 = cn n

[1 − (−1) ]

Z1

Pn (cos θ) d cos θ = cn

cos α

2 2n + 1 2 2n + 1

(4.38)

Het enige probleempunt ter bepaling van de coëfficiënten cn is de berekening van de integraal in vgl. (4.38). We gebruiken hiervoor de volgende recursiebetrekking [47]: (2n + 1) Pn (x) =

d [Pn+1 (x) − Pn−1 (x)] dx

(4.39)

Algemeen geldt dus dat: Zx2 Pn (x) dx =

x1

1 [(Pn+1 (x2 ) − Pn−1 (x2 )) − (Pn+1 (x1 ) − Pn−1 (x1 ))] 2n + 1

(4.40)

Toepassing hiervan in vgl. (4.38) levert: cn

= (4.25)

=

2n + 1 1 − (−1)n [(Pn+1 (1) − Pn−1 (1)) − (Pn+1 (cos α) − Pn−1 (cos α))] 2n + 1 2 1 − (−1)n [Pn−1 (cos α) − Pn+1 (cos α)] (4.41) 2

Gebruik makend van de coëfficiënten van de Fourier-Legendre-reeks kunnen we de onbekende an bepalen dankzij de randvoorwaarde zoals beschreven in vgl. (4.31): J (cos θ)

(4.31)

=

σ

∞ X

n an Rn−1 Pn (cos θ)

n=0 (4.41)

=

∞ X 1 − (−1)n

n=0

2

[Pn−1 (cos α) − Pn+1 (cos α)] Pn (cos θ)

(4.42)

Wegens de orthogonaliteit moet deze laatste betrekking opnieuw gelden voor elke afzonderlijke term: σ n an Rn−1 = an =

1 − (−1)n [Pn−1 (cos α) − Pn+1 (cos α)] 2 1 − (−1)n 1 [Pn−1 (cos α) − Pn+1 (cos α)] 2 σ n Rn−1

∀n ∈ N0 (4.43)

55

Hoofdstuk 4. Numerieke analyse van numerieke methoden De oplossing van deze analytische randsituatie wordt dus: ∞ X 1 − (−1)n

1 V (r, θ) = [Pn−1 (cos α) − Pn+1 (cos α)] rn Pn (cos θ) 2 σ n Rn−1 n=1 ∞ X 1 1 [P2m (cos α) − P2m+2 (cos α)] r2m+1 P2m+1 (cos θ) = 2m + 1 σ R2m

(4.44) (4.45)

m=0

Het werd doorheen de afleiding duidelijk dat deze oplossing een veralgemening is van het geval met een delta-bron (vgl. (4.30)). Dit geval is terug te vinden indien we cos α → 1 in fig. 4.2. We bekomen in de limiet een delta-functie in plaats van het rechthoekig profiel. Berekenen we: lim [P2m (cos α) − P2m+2 (cos α)]

cos α→1

= (4.25)

=

=

lim [(1 − P2m+2 (cos α)) − (1 − P2m (cos α))]

cos α→1

lim

Z1

0 0 P2m+2 (cos θ) − P2m (cos θ) d cos θ

cos α→1 cos α 0 P2m+2 (1) −

0 P2m (1)

(4.46)

Gebruiken we opnieuw de recursiebetrekking zoals in (4.39) dan is: 0 0 P2m+2 (x) − P2m (x) = (4m + 3) P2m+1 (x)

∀x ∈ [−1, 1]

(4.47)

Hierdoor bedraagt: lim [P2m (cos α) − P2m+2 (cos α)]

cos α→1

(4.46)

=

(4.47)

=

(4.25)

=

0 0 P2m+2 (1) − P2m (1)

(4m + 3) P2m+1 (1) 4m + 3

(4.48)

Substitutie hiervan in de oplossing (4.45) van het bolkapprobleem levert inderdaad de oplossing met de delta-bron op ((4.30)).

4.2.2

Numerieke simulatie

Om de potentiaal numeriek te bepalen hanteren we de numerieke methoden zoals beschreven in sectie 3.3.2. We gebruiken een kubisch rooster bestaande uit 101 × 101 × 101 voxels met ∆x = ∆y = ∆z = 0.001 m. Rond het centrum met roostercoördinaten (51, 51, 51) plaatsen we een sfeer met een straal van 30 voxels. Deze sfeer heeft een elektrische geleidbaarheid van S 0.5 m . Zoals aangeduid in fig. 4.1 wordt elke naald in één voxel ingebracht. De naalden zelf worden langs de y-richting ingebracht. De bekomen numerieke waarden voor de potentiaal worden weergegeven in fig. 4.3. In het linkerbovenluik duiden we de convergentieratio aan ten opzichte van het aantal iteraties nodig in de Preconditioned BiCGSTAB-routine. We zien dat na 350 iteraties een convergentieratio van 10−16 wordt bereikt. Dit zou normaliter tot nauwkeurige resultaten moeten leiden. In de overige deelfiguren van fig. 4.3 wordt telkens de potentiaal weergegeven in de snede

56

0

z [voxelnummers]

log(convergentieratio)


−5 −10 −15 −20 0

100

200

300

100 80 60

0 40 20 20

400

80

−0.005

60

−0.01

40

−0.015

20

−0.02 20

40

60

80

100

x [voxelnummers]

y [voxelnummers]

z [voxelnummers]

0

y = 51

40

60

80

100

−50

y [voxelnummers]

Aantal iteraties [-] 100

50

x = 51

100

50

z = 51

80 60 0 40 20 20

40

60

80

100

−50

x [voxelnummers]

Figuur 4.3: Numerieke resultaten van een elektromagnetisch deelprobleem.

x = 51, y = 51 of z = 51 (uitgedrukt in roostercoördinaten). In de x- en z-snede valt onmiddellijk op te merken dat de potentiaal heel gelokaliseerd is. De regio’s met hoge potentialen zijn gesitueerd in de nabijheid van de voxel waar de naalden ingebracht zijn. In de y-snede zien we dat de potentiaal quasi overal 0 is. Dit valt te verklaren op basis van de symmetrie van de simulatieomgeving. De ene naald injecteert stroom, terwijl de andere naald stroom extraheert. Deze punten liggen symmetrisch ten opzichte van de beschouwde snede y = 51. De ene naald ageert dus als een spiegelbron voor de andere naald. De potentiaal moet dan identisch nul zijn in de volledige snede. In het linkerbenedenluik van fig. 4.3 zijn er desalniettemin kleine afwijkingen door o.a. discretisatiefouten. We mogen besluiten dat onze numerieke methode dus in staat is om resultaten te genereren die voldoen aan de symmetrie van het probleem. 4.2.2.1

Overeenstemming met delta-bron

We kunnen nu deze numerieke resultaten vergelijken met onze analytische formule afgeleid voor een delta-bron. Concreet vergelijken we de genormeerde potentialen bekomen met vgl. (4.30) met de genormeerde potentialen van de numerieke simulatie en dit in een snede z = 51. Beide resultaten worden weergegeven in fig. 4.4. Op basis van deze figuur kunnen we opnieuw besluiten dat het potentiaalverloop heel sterk gelokaliseerd is ter hoogte van de naaldposities. Desalniettemin merken we op dat beide oplossingen niet volledig overeenstemmen. in fig. 4.5. Hier wordt de Dit wordt bevestigd Vnumeriek −Vanalytisch logaritme van het relatieve verschil log uitgezet in de snede z = 51. Vanalytisch In de nabijheid van de naalden vinden we vooral een gele regio terug die overeenstemt met een relatieve afwijking van circa 10%. In het overgrote merendeel van de snede moeten we echter vaststellen dat er een relatieve afwijking van 10−0.5 of dus 30% is. Dit valt te verklaren door het feit dat een delta-bron slechts stroom injecteert in een infinitesimaal punt, terwijl de numerieke simulatie stroom injecteert in een voxel met een eindig volume. Onze analytische uitdrukking is dus logischerwijze een vrij ruwe benadering van de numerieke simulatie.

57


Figuur 4.4: Vergelijking van de genormeerde numerieke en analytische oplossingen in het geval van de delta-bron.

100

0

90

y [voxelnummers]

80

−0.5

70 −1

60 50

−1.5

40 30

−2

20 10

−2.5 20

40

60

80

100

x [voxelnummers] Figuur 4.5: Logaritme van het relatieve verschil tussen de genormeerde numerieke en analytische oplossingen in het geval van de delta-bron.

4.2.2.2

Overeenstemming met bolkapelektrode

Om enige ruimtelijke uitgestrektheid te verkrijgen, ontwikkelden we het analytische model van de bolkapelektrode. Zoals in de vorige subsectie zullen we de genormeerde analytische potentialen uit vgl. (4.45) vergelijken met de potentialen verkregen uit de numerieke simulatie. Beide resultaten in de snede z = 51 worden weergegeven in fig. 4.6. Opnieuw moeten we opmerken dat de oplossingen niet volledig overeenstemmen. De relatieve verschillen zijn echter kleiner geworden zoals afgebeeld in fig. 4.7. In de buurt van de naalden bedragen de relatieve verschillen tussen 10−2.5 en 10−3.5 . Procentueel betekent dit een afwijking van slechts 0.3 − 3%. Het overige deel van de snede vertoont een relatieve afwijking tussen 10 en 17%. Ten opzichte van de situatie met de delta-bron, betekent dit dat onze analytische oplossing in betere mate overeenstemt met de numerieke

58


Figuur 4.6: Vergelijking van de genormeerde numerieke en analytische oplossingen in het geval van de bolkapelektrode.

100 −0.5

y [voxelnummers]

80 −1 60 −1.5

−2

40

−2.5 20 −3 20

40

60

80

100

x [voxelnummers]

Figuur 4.7: Logaritme van het relatieve verschil tussen de genormeerde numerieke en analytische oplossingen in het geval van de bolkapelektrode.

resultaten. Een volledige overeenstemming wordt echter niet teruggevonden. Dit was te verwachten aangezien de bolkapelektrode geen ruimtelijke uitgestrektheid heeft in de radiale richting. Inderdaad, de bolkapelektrode injecteert stroom ter hoogte van het oppervlakte r = R, terwijl de numerieke simulatie stroom injecteert in één volledige voxel waardoor er wel degelijk een uitgestrektheid is in de radiale richting. De bolkapelektrode vormt dus weerom een benadering van de numerieke situatie. We mogen echter besluiten dat de potentialen bekomen in de numerieke simulatie correct zijn. Vooreerst is onze numerieke methode in staat om de symmetrievoorwaarden van het elektromagnetische probleem correct te reproduceren. Verder is er een goede overeenstemming met de analytische resultaten in het geval van de bolkapelektrode ondanks het feit dat deze slechts een benadering vormt van de numerieke simulatie.


4.3 4.3.1

59

Thermisch deelprobleem Analytische limietsituatie

In deze sectie trachten we onze numerieke methode aangaande het thermische deelprobleem te analyseren. We vertrekken hierbij van het wiskundige model dat dit deelprobleem beschrijft. We hernemen vgl. (3.26): ρ (~r) c (~r)

∂T (~r, t) ∂t

= ∇ · [k (~r, T (~r, t)) ∇T (~r, t)] − ωρb cb (T (~r, t) − Tc ) 2 ~ + σ (~r, T (~r, t)) E (~r, t)

(4.49)

De timestepping via de modified midpoint-methode kan geverifieerd worden door de analytische oplossing te beschouwen bij het aanleggen van een eenheidspuls op t = 0 en ~x = ~0 (T (x, y, z, t = 0) = δ (x) δ (y) δ (z)) in vgl. (4.49)). Deze analytische oplossing kan men dan vergelijken met de numerieke resultaten in dit limietgeval. Aangezien de correctheid van het elektrisch veld in de vorige sectie gevalideerd wordt, mogen we ons beperken tot het wis~ ~ kundige model zonder de term die het elektrisch vermogen beschrijft E = 0 . Een verdere niet-essentiële vereenvoudiging bestaat erin Tc = 0 te stellen. Ten slotte werken we in een uniform medium zodat ∇k = ~0. Zo verkrijgen we het volgende beginwaardenprobleem:  ρc ∂T (x,y,z,t) = k∇2 T (x, y, z, t) − ωρ c T (x, y, z, t) b b ∂t (4.50) T (x, y, z, t = 0) = δ (x) δ (y) δ (z)

Triviale uitwerking levert:   ∂T (x,y,z,t) − k ∇2 T (x, y, z, t) + ωρb cb T (x, y, z, t) = 0 ∂t ρc ρc T (x, y, z, t = 0) = δ (x) δ (y) δ (z)

Stellen we K = als:

k ρc

en A =

ωρb cb ρc ,

(4.51)

dan is het uiteindelijke beginwaardenprobleem te schrijven

  ∂T (x,y,z,t) − K∇2 T (x, y, z, t) + AT (x, y, z, t) = 0 ∂t T (x, y, z, t = 0) = δ (x) δ (y) δ (z)

(4.52)

Beschouwen we nu de volgende transformatie:

T (x, y, z, t) = T 0 (x, y, z, t) exp (−At)

(4.53)

dan transformeert beginwaardenprobleem (4.52) in het volgende beginwaardenprobleem:  ∂T 0 (x, y, z, t)   exp (−At) − AT 0 (x, y, z, t) exp (−At) − K exp (−At) ∇2 T 0 (x, y, z, t)   ∂t + AT 0 (x, y, z, t) exp (−At) = 0     0 T (x, y, z, t = 0) = δ (x) δ (y) δ (z) (4.54) Dit resulteert in de omschrijving van de fundamentele oplossing van de warmtevergelijking in een homogeen medium:   ∂T 0 (x,y,z,t) − K∇2 T 0 (x, y, z, t) = 0 ∂t (4.55) T 0 (x, y, z, t = 0) = δ (x) δ (y) δ (z)

60


De oplossing van dit gekende beginwaardenprobleem is (zie sectie A.2 voor de afleiding van dit resultaat): 2 1 x + y2 + z2 0 T (x, y, z, t) = − H (t) (4.56) 3 exp 4Kt (4πKt) 2 Via transformatie (4.53) verkrijgen we de volgende uitdrukking voor het temperatuursprofiel: 2 1 x + y2 + z2 T (x, y, z, t) = − exp (−At) H (t) (4.57) 3 exp 4Kt (4πKt) 2 Het temperatuursprofiel uitgedrukt in de oorspronkelijke variabelen is dan: ρc 3 x2 + y 2 + z 2 ωρb cb 2 exp −ρc exp − T (x, y, z, t) = t H (t) 4πkt 4kt ρc

(4.58)

Zoals in de numerieke analyse van het elektromagnetisch deelprobleem aangekaart, is een fundamentele oplossing niet zo handig om te vergelijken met een numerieke simulatie. In plaats van de bronterm δ (x) δ (y) δ (z) in vgl. (4.50) leggen we een blokfunctie aan als bron: u (x) u (y) u (z). We definiëren u (xi ) hierbij als:  1 if |x | ≤ n ∆xi i 2 (4.59) u (xi ) = 0 if |x | > n ∆xi i

2

Voor de duidelijkheid geven we het functievoorschrift van u (xi ) ook visueel weer in fig. 4.8.

(xi) 1

-n

Δxi 2

n

Δxi 2

xi

Figuur 4.8: Grafische weergave van het functievoorschrift voor u (xi ).

In distributionele zin bekomen we dus het volgende beginwaardenprobleem: ∂T (x, y, z, t) − K∇2 T (x, y, z, t) + AT (x, y, z, t) = u (x) u (y) u (z) δ (t) ∂t

(4.60)

Gelukkig kunnen we nog steeds gebruik maken van de fundamentele oplossing om vgl. (4.60) op te lossen. De oplossing van dit probleem is immers niets anders dan de convolutie van de

61


fundamentele oplossing met u (x) u (y) u (z). We bewijzen deze uitspraak als volgt. Beschouw de lineaire differentiaaloperator L. Dan voldoet de fundamentele oplossing F aan LF (x) = δ (x). Dan is: L (F ∗ G) = (LF (x)) ∗ G = δ (x) ∗ G = G (x)

(4.61)

of dus is de convolutie van de fundamentele oplossing met een andere functie de oplossing van het beginwaardenprobleem met die tweede functie als bronterm. Vooreerst herschrijven we de fundamentele oplossing (4.57) als: " 3 # Y 0 T1 (x, y, z, t) = T1 (xi , t) exp (−At) H (t) (4.62) i=1

met:

x2i = − 1 exp 4Kt (4πKt) 2 Dan kunnen we de oplossing van beginwaardenprobleem (4.60) berekenen als volgt: ZZZ T (x, y, z, t) = d~x0 T1 ~x − ~x0 , t u ~x0 1

T10 (xi , t)

(4.63)

R3

(4.62)

=

+∞ +∞ +∞ Z Z Z 0 0 dz 0 T10 x − x0 , t T10 y − y 0 , t T10 z − z 0 , t dy dx

(4.59)

=

−∞

−∞ 0

−∞

u x u y u z 0 exp (−At) H (t)    n ∆x n ∆y Z2 Z2     T10 y − y 0 , t dy 0  T10 x − x0 , t dx0   ·

   

0

−n ∆x 2



 ·  =

n ∆z 2

Z

−n ∆z 2



n

−n ∆y 2



 T10 z − z 0 , t dz 0   exp (−At) H (t)

∆xi 2

Z 3 Y  

i=1 ∆x −n 2 i



 T10 xi − x0i , t dx0i   exp (−At) H (t)

(4.64)

De enige stap die nog nodig is om tot de oplossing te komen, is dus de berekening van de volgende integraal: n

∆xi

Z2

(4.63) xi − x0i , t dx0i =

T10

∆x −n 2 i

n

∆xi

Z2

∆x −n 2 i

"

(xi − x0i )2 1 exp − 4Kt (4πKt) 2 1

#

dx0i

√ x0 −x i Stellen we u gelijk aan √i4Kti , dan is dx0i = 4Kt du. Als x0i = ±n ∆x 2 dan is u = Substitutie hiervan in uitdrukking (4.65) levert: n

n

∆xi 2

Z

∆x −n 2 i

T10

xi −

x0i , t

dx0i

1 =√ π

∆xi −xi √2 4Kt

Z

∆x −n 2 i −xi √ 4Kt

exp −u2 du

(4.65)

∆x

±n 2 i −xi √ . 4Kt

(4.66)

62


Met behulp van de errorfunctie kunnen we dit resultaat herschrijven. De errorfunctie is een speciale wiskundige functie die gedefinieerd is als: 2 erf (x) = √ π

Zx 0

exp −t2 dt

(4.67)

Gebruiken we deze functie in vgl. (4.66), dan verkrijgen we: n

∆xi

Z2

−n

T10

xi −

∆xi 2

x0i , t

dx0i

Aangezien erf (−x) = −erf (x) n

∆xi

Z2

−n

∆xi 2

T10

" 1 = erf 2

i n ∆x − xi √2 4Kt

!

− erf

i − xi −n ∆x √2 4Kt

!#

(4.68)

∀x ∈ R, bekomen we:

" 1 erf xi − x0i , t dx0i = 2

i xi + n ∆x 2 √ 4Kt

!

− erf

i xi − n ∆x 2 √ 4Kt

!#

(4.69)

Vervangen we dit deelresultaat in vgl. (4.64), dan kunnen we de oplossing van beginwaardenprobleem (4.60) schrijven als: ! !## " " 3 i i Y xi − n ∆x xi + n ∆x 1 2 2 √ √ − erf exp (−At) H (t) (4.70) erf T (x, y, z, t) = 2 4Kt 4Kt i=1

4.3.2

Numerieke simulatie

Beschouwen we nu een homogene omgeving van gezond weefsel wiens thermische eigenschappen samengevat worden in tabel 3.1. We delen deze omgeving op in 51 × 51 × 51 kubische voxels waarvoor ∆x = ∆y = ∆z = 0.001 m. Het centrum van dit kubisch rooster is dus gevestigd in het roosterpunt met roostercoördinaten (26, 26, 26). We stellen nu de begintemperatuur van de voxels met roostercoördinaten (i, j, l) die voldoen aan 16 ≤ i ≤ 36, 16 ≤ j ≤ 36, 16 ≤ l ≤ 36 op 100 K. De overige voxels buiten dit kubisch deelgebied hebben een begintemperatuur van 0 K. We wensen in deze omgeving de analytische oplossing bepalen van de volgende partiële differentiaalvergelijking: ∂T (x, y, z, t) − K∇2 T (x, y, z, t) + AT (x, y, z, t) = 0 ∂t

(4.71)

met de hierboven besproken beginvoorwaarden. Deze situatie vertoont een grote analogie met de analytische situatie van de vorige sectie. Inderdaad de begintemperaturen voldoen aan: T (x, y, z, 0) = 100 u (x) u (y) u (z) (4.72) met n = 21 in de definitie van u (xi ) (vgl. (4.59)). De analytische oplossing van dit probleem volgt dan ook direct uit de oplossing beschreven in vgl. (4.70): " " ! !## 3 ∆xi i Y xi + 21 ∆x x − 21 1 i 2 2 √ √ T (x, y, z, t) = 100 erf − erf exp (−At) H (t) (4.73) 2 4Kt 4Kt i=1

63


Nu kunnen we deze analytische oplossing vergelijken met de numerieke resultaten voor deze realistische simulatieomgeving. Op basis hiervan kunnen we een uitspraak doen over de nauwkeurigheid van onze implementatie. Vooreerst bekijken we de temperatuur in één voxel gedurende een tijdsverloop van 300 s. In fig. 4.9 wordt de centrale voxel (26, 26, 26) beschouwd. In deze figuur is duidelijk te merken 100

Numerieke simulatie Analytische oplossing

90

Temperatuur [K]

80 70 60 50 40 30 20 10 0 0

50

100

150

200

250

300

Tijd [s]

Figuur 4.9: Numerieke en analytische oplossing in het punt (26, 26, 26) binnen de centrale regio.

dat de analytische en numerieke oplossing elkaar heel goed benaderen. Dit wordt verder bevestigd in fig. 4.10. Hierin wordt de verhouding van de numerieke en analytische oplossing in het centrum weergegeven.

Temperatuursverhouding [-]

1.003

1.002

1.001

1

0.999

0.998

0.997 0

50

100

150

200

250

300

Tijd [s]

Figuur 4.10: Verhouding van de numerieke en analytische oplossing in het punt (26, 26, 26) binnen de centrale regio.

We zien dat in de centrale regio de afwijking van het numerieke resultaat ten opzichte van de analytische oplossing nauwelijks 0.2% bedraagt. Het is echter ook belangrijk om eens te

64


kijken naar een voxel in het gebied met begintemperatuur 0 K. Dit gebied is veel gevoeliger aan numerieke afwijkingen door bijvoorbeeld discretisatiefouten. Het temperatuursverloop van een voxel gelokaliseerd ter hoogte van (11, 11, 11) wordt afgebeeld in fig. 4.11. 0.7

Numerieke simulatie Analytische oplossing

Temperatuur [K]

0.6 0.5 0.4 0.3 0.2 0.1 0 0

50

100

150

200

250

300

Tijd [s]

Figuur 4.11: Numerieke en analytische oplossing in het punt (11, 11, 11) buiten de centrale regio.

We zien dat de numerieke oplossing de analytische oplossing benadert, maar her en der moeten we een kleine afwijking tolereren. Het is logisch dat de abrupte overgang van 100 K naar 0 K een grotere fout zal veroorzaken dan in de centrale regio waar de overgang veel zachter is tussen aanleunende voxels. 1.4


1.2 1 0.8 0.6 0.4 0.2 0 0

50

100

150

200

250

300

Tijd [s]

Figuur 4.12: Verhouding van de numerieke en analytische oplossing in het punt (11, 11, 11) buiten de centrale regio.

Kijken we naar de temperatuursverhouding (fig. 4.12), dan zien we dat er een heel slechte overeenkomst is tot circa 16 s. Dit wordt veroorzaakt door het feit dat de analytische oplossing in dit gebied tussen de 10−17 en 10−5 K bedraagt. Een kleine absolute afwijking door afrondfouten zal aldus een grote relatieve fout veroorzaken. Deze grote relatieve afrondfouten

65


hebben echter geen invloed op de latere tijdstippen. In het tijdsgebied tussen 50 s en 300 s is de maximale afwijking 4.6%. Vervolgens bestuderen we het temperatuursverloop in een spatiale doorsnede op één enkel tijdstip. We kiezen ervoor om op t = 100 s de snede {y = 0, z = 0} te bekijken. In roosterco¨ ordinaten komt deze snede overeen met j = 26 en l = 26 en waarbij i kan variëren tussen 1 en 51 of equivalent x ∈ [−0.025, 0.025] m. Het temperatuursprofiel in de desbetreffende snede wordt weergegeven in fig. 4.13. 45 Numerieke simulatie Analytische oplossing

40

Temperatuur [K]

35 30 25 20 15 10 5 0 −0.025

−0.02

−0.015

−0.01

−0.005

0

0.005

0.01

0.015

0.02

0.025

x [m]

Figuur 4.13: Numerieke en analytische oplossing op het tijdstip t = 100 s in de snede {y = 0, z = 0}.

Het is opnieuw duidelijk dat de numerieke simulatie in staat is om correcte resultaten te genereren. Een blik op de temperatuursverhouding (fig. 4.14) leert ons dat deze indruk juist is. 1.05


1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 −0.025

−0.02

−0.015

−0.01

−0.005

0

0.005

0.01

0.015

0.02

0.025

x [m]

Figuur 4.14: Verhouding van de numerieke en analytische oplossing op het tijdstip t = 100 s in de snede {y = 0, z = 0}.

Op de uiteinden van het beschouwde domein blijkt de relatieve fout opnieuw vrij groot te zijn.

66


De reden is dezelfde als in het temporele geval: de absolute waarde van de oplossingen op die posities is dermate klein dat een kleine absolute fout grote relatieve fouten kan genereren. In het tussengebied x ∈ [−0.020, 0.020] m bedraagt de maximale relatieve fout nauwelijks 0.15%. We kunnen dus zowel temporeel als spatiaal gewag maken van een nauwkeurige numerieke methode aangaande het thermisch deelprobleem. Ten slotte kunnen we het stabiliteitsgedrag van onze numerieke methode controleren. In vgl. (3.51) hebben we een bovengrens voor de tijdstap ∆t afgeleid. Substitueren we de materiaaleigenschappen van gezond weefsel (tabel 3.1) en de fysische constanten van de bloedperfusieterm (tabel 3.2) in deze uitdrukking, dan verkrijgen we: ∆t ≤

4 · 0.53 ·

2 · 1016 · 3500 = 1.113555 s + 0.0064 · 1000 · 4180

(4.74)

3 0.0012

Indien we een tijdstap ∆t kiezen die kleiner is dan deze bovengrens zou onze numerieke methode, althans in theorie, stabiel moeten zijn. We simuleren hetzelfde probleem als hierboven met verschillende tijdstappen en analyseren de bekomen resultaten. Zo kunnen we het temperatuursverloop in de centrale voxel bekijken voor verschillende tijdstappen. Dit wordt in fig. 4.15 weergegeven. 100 Analytische oplossing ∆ t = 1.3043 s ∆ t = 1.333 s ∆ t = 1.3636 s ∆ t = 1.5 s ∆t=2s

90

Temperatuur [K]

80 70 60 50 40 30 20 10 0 0

50

100

150

200

250

300

Tijd [s]

Figuur 4.15: Numerieke en analytische oplossing in het punt (26, 26, 26) voor verschillende stapgrootten in de tijd.

Het is onmiddellijk duidelijk dat een ∆t = 1.3043 s (230 tussenstappen voor een simulatietijd van 300 s), opnieuw heel dicht in de buurt zit van de analytische oplossing. Met deze stapgrootte in de tijd verkrijgen we met andere woorden een stabiele methode. Nemen we een grotere stapgrootte, bv. ∆t = 1.333 s, dan divergeert de numerieke oplossing van de analytische oplossing na ongeveer 170 s. De absolute fout neemt exponentieel toe zoals men kan opmerken uit de steile toename van de divergentie. In fig. 4.15 is duidelijk merkbaar dat voor nog grotere stapgroottes eveneens instabiliteiten optreden. Hoe groter de tijdstap, des te sneller de afwijking tevoorschijn zal komen. We besluiten dat een ∆t = 1.3043 s de limiet voor stabiliteit is. Deze bovengrens is groter dan de bovengrens die we theoretisch voorspeld hebben in vgl. (4.74). Dit valt te verwachten

67


aangezien de theoretische bovengrens is afgeleid voor de FTCS Euler-methode. De modified midpoint-methode is blijkbaar iets stabieler dan deze standaardmethode. Desalniettemin kunnen we concluderen dat als we voldoen aan de theoretische bovengrens wat betreft stapgrootte in de tijd, we een stabiele numerieke methode verkrijgen. Dit is een heel belangrijk tussenresultaat, aangezien we een voorwaarde verkregen hebben die ons toelaat om de numerieke berekeningen op een correcte manier uit te voeren. Om terug te komen op de instabiliteiten bij grotere tijdstappen geven we de temperatuur in een dwarsdoorsnede weer op verschillende tijdstippen (fig. 4.16). −9

−8

x 10

50

40

8

40

30

6

20

4

10

2

20

30

40

50

x [voxelnummers]

12 10

30

8 6

20

4 10

2

0 10 x 10

50

2.5

t = 19.5 s 40

2

30

1.5

20

1

10

0.5

20

30

40

50

x [voxelnummers]

−6

50

14

t = 15 s y [voxelnummers]

t = 12 s

10

y [voxelnummers]

x 10

10

0 34

x 10

t = 300 s

3 2

40

y [voxelnummers]

y [voxelnummers]

50

1 30 0 20

−1 −2

10

−3 0 10

20

30

x [voxelnummers]

40

50

10

20

30

40

50

x [voxelnummers]

Figuur 4.16: Numerieke oplossing in de snede z = −0.025 m waarbij de evolutie van de numerieke instabiliteit wordt gevisualiseerd.

Indien we ∆t = 1.5 s gebruiken, zien we correcte resultaten tot circa t = 15 s. Vanaf dan ontstaat er in de vier hoekpunten van de rode zone een dambordpatroon van kleine fouten (zie bv. t = 19.5 s). De absolute waarde van deze fouten is betrekkelijk klein (orde 10−6 ), wat ook gereflecteerd wordt in het feit dat de numerieke oplossing nog steeds de analytische oplossing benadert. Deze fouten zullen echter exponentieel toenemen en door het convectiemechanisme verspreid worden over de hele simulatieomgeving. Op het einde van de simulatie (t = 300 s) is de numerieke oplossing compleet onrealistisch, daar ze over de hele snede temperaturen aangeeft tussen −3 · 1034 en 3 · 1034 K. Dit onderstreept nogmaals het belang van een theoretische bovengrens om dergelijke onrealistische situaties te vermijden. Indien we geen aandacht zouden verlenen aan het onderwerp van de stabiliteit, zouden we verkeerde conclusies kunnen trekken uit de bekomen numerieke resultaten.

Hoofdstuk 5

Simulatie van thermotherapie However beautiful the strategy, you should occasionally look at the results. Winston Churchill

5.1

Configuratie van de simulatie

Onder de configuratie van de simulatie verstaan we de geometrie, de materiaalparameters en de beslissingscriteria die gebruikt worden voor de simulaties in dit hoofdstuk. Bij het voorstellen van de geometrie gebruiken we opnieuw kubusvormige voxels met afmetingen ∆x = ∆y = ∆z = 0.001 m. We opteren ervoor om een grote omgeving te simuleren met behulp van 101 × 101 × 101 voxels. In fig. 5.1 wordt onze simulatieomgeving visueel weergegeven door middel van een centrale dwarsdoorsnede. 100 90

y [voxelnummers]

80 70 60 50 40 30 20 10 20

40

60

80

100

x [voxelnummers]

Figuur 5.1: Dwarsdoorsnede van de simulatieomgeving.

Deze figuur geeft een dwarsdoorsnede weer van een sferisch object (de lever) met materiaaleigenschappen die overeenstemmen met goedaardig weefsel. De straal van deze bol bedraagt 40 voxels. In deze bol bevindt zich echter ook een sferisch volume tumorweefsel met een 68

69

Hoofdstuk 5. Simulatie van thermotherapie

straal van 15 voxels. De naalden bevinden zich in de tumor om daar het weefsel lokaal te verwarmen. In de tip van de bovenste naald zal er stroom ge¨ınjecteerd worden, terwijl er in de tip van de onderste naald extractie van stroom zal plaatsvinden. Naast een dwarsdoorsnede kunnen we in MATLAB ook een impressie geven van hoe deze simulatieomgeving eruit ziet in drie dimensies. Een combinatie van verschillende sneden wordt afgebeeld in fig. 5.2.

z [voxelnummers]

100 80 60 40 20

100 20

80 40

x [voxelnummers]

60 60

40 80

y [voxelnummers]

20 100

Figuur 5.2: Driedimensionale impressie van de simulatieomgeving.

Verder geven we de materiaaleigenschappen van de verschillende weefsels weer in tabel 5.1. Tabel 5.1: Materiaaleigenschappen van de verschillende weefsels die aanwezig zijn in de simulatieomgeving.

Type weefsel

lucht

gezond weefsel

maligne weefsel

Voxellabel σ0 S m−1 σ1 K −1 k0 W m−1 K −1 k1 W m−1 K −2 c J kg −1 K −1 ρ kg m−3

1 10−7 0 0.025 0 3500 1016

2 0.36 0.015 0.53 0.00116 3500 1016

3 0.45 0.015 0.53 0.00116 3500 1016

Het grote verschil ten opzichte van de waarden in tabel 3.1 is de veranderde ρ en c voor lucht. Voor lucht gebruiken we dezelfde massadichtheid en specifieke warmtecapaciteit als die voor gezond en maligne weefsel. We voeren met andere woorden een soort homogenisatie uit. Op het eerste zicht lijkt een dergelijke propositie te gek voor woorden. Dit is immers een grove benadering van de fysische werkelijkheid. Het grote voordeel van deze benadering ligt hem in het feit dat onze bovengrens voor de stabiele tijdstap van een aanvaardbare grootte-orde blijft. Zouden we de werkelijke mate-


70

riaaleigenschappen voor lucht gebruiken, dan moeten we krachtens vgl. (3.51) een tijdstap nemen kleiner dan: ∆t ≤

4 · 0.025 ·

2 · 1 · 1012 = 0.00674667 s + 0.0064 · 1000 · 4180

3 0.0012

(5.1)

Een dergelijke tijdstap is allesbehalve interessant aangezien we per seconde dat we RFA toepassen dan 149 tijdstappen moeten gebruiken. Dit resulteert in een computationeel totaal inefficiënte methode. Bovendien tonen simulaties met deze benadering aan dat de numerieke waarden voor de temperatuur in de gezonde en maligne weefsels niet verschillen van de werkelijke waarden. Dit hoeft ons niet te verwonderen aangezien lucht een heel kleine thermische geleidbaarheid heeft en de warmtebronnen gelokaliseerd zijn in het tumorweefsel. Warmte ontstaat binnen de tumoren en propageert naar buiten toe. De invloed van lucht heeft dus slechts een minuscuul effect op de temperatuur in de weefsels waar onze interesse naar uitgaat. Deze invloed is zodanig klein dat het niet observeerbaar is omdat afrondfouten sowieso de resolutie van de numerieke waarden beperkt. Een belangrijk aspect van een simulatie van RFA is het beslissingscriterium omtrent de stroominjectie. Het lijkt ons niet opportuun om gedurende de hele RFA continu stroom te injecteren aangezien de maximumtemperatuur in het weefsel tot temperaturen hoger dan de maximumtemperatuur van circa 100 ◦ C zou kunnen oplopen. In de literatuur beschrijft men vaak een procedure voor een monopolaire configuratie. Deze verschilt van auteur tot auteur, maar meestal opteert men ervoor om de temperatuur te laten oplopen tot 100 ◦ C en dan de temperatuur te laten zakken tot een niet nader bepaalde temperatuur waarop er geen schade wordt aangericht. Wij kiezen er in dit hoofdstuk voor om eveneens een temperatuurscriterium te hanteren. We injecteren een vaste stroom van 0.07 A tot de maximumtemperatuur in het weefsel 80 ◦ C bedraagt. Daarna zetten we de stroominjectie stop tot de maximumtemperatuur gedaald is tot 50 ◦ C. Op dat moment schakelen we de stroom weer aan zodat de hele cyclus opnieuw begint. De reden waarom we 80 ◦ C verkiezen boven 100 ◦ C ligt hem in fig. 2.5. Deze figuur toont aan dat er irreversibele schade optreedt boven een temperatuur van 80 ◦ C. Laten we de temperatuur dan dalen, dan kunnen we de mathematische uitdrukkingen voor de materiaalparameters niet langer gebruiken. Het is bovendien onmogelijk om de materiaaleigenschappen te beschrijven van weefsels die herhaaldelijk irreversibele schade ondergaan hebben. Dit schijnt de meeste auteurs in de literatuur niet te deren, ondanks het feit dat de correctheid van hun simulaties hierdoor op losse schroeven komt te staan. Verder motiveren we de ondergrens van 50 ◦ C door te verwijzen naar subsectie 2.1.3 waarin deze temperatuur als ondergrens wordt beschouwd van het temperatuursgebied waar we thermische schade kunnen aanrichten.

5.2

Invloed van thermische en elektrische parameters

Eens de simulatieomgeving is vastgelegd, gaan we na wat de invloed van de verschillende weefsels en hun specifieke thermische en elektrische eigenschappen is op de temperatuursdistributie.

71


Vooreerst kunnen we bekijken wat er verandert aan de temperatuur indien we het tumorweefsel verwijderen uit de simulatieomgeving van de vorige sectie. We simuleren zowel de configuratie met een tumor, als een configuratie zonder tumor gedurende één cyclus of puls (1 maal stroom aanschakelen en afschakelen). De temperatuur op de positie (51, 45, 51) nabij de naalden wordt weergegeven in fig. 5.3 voor het volledige tijdsverloop. 355

zonder tumor met tumor

350

Temperatuur [K]

345 340 335 330 325 320 315 310 0

1

2

3

4

5

6

Tijd [s]

Figuur 5.3: Vergelijking van het temperatuursverloop in het punt (51, 45, 51) in aan- of afwezigheid van een tumor.

Het is onmiddellijk duidelijk dat er vrij grote verschillen zijn tussen de twee configuraties. Door de verschillende elektrische geleidbaarheden van goedaardig en maligne weefsel zullen de potentialen andere waarden aannemen. De logaritme van het verschil in potentiaal (log |∆φ|) wordt weergegeven in fig. 5.4. 100 0.5

90

0

z [voxelnummers]

80

−0.5

70

−1

60

−1.5 50 −2 40 −2.5 30 −3 20 −3.5 10

−4 20

40

60

80

100

y [voxelnummers]

Figuur 5.4: Logaritme van het verschil in potentiaal tussen beide configuraties in de snede x = 51.

Het grootste verschil in potentiaal bevindt zich in de tumorregio. Desalniettemin zijn er

72


eveneens relatief grote verschillen in het gezonde weefsel. Blijkbaar beperkt de invloed van het tumorweefsel zich niet alleen tot de regio waar de tumor resideert. Dit verschil in potentiaal heeft een grote invloed op de elektrische bronterm in de Pennes bio-warmtevergelijking. Hierdoor zullen de temperaturen verschillen in beide configuraties. Bovendien zal hierdoor de thermische geleidbaarheid k (T ) verschillen wat eveneens een belangrijk effect heeft op de Pennes bio-warmtevergelijking. Doordat de eerste configuratie (zonder tumor) sneller de maximumtemperatuur bereikt, zal de stroom sneller afgeschakeld worden. Dit resulteert dan in de dalende temperatuur zoals afgebeeld in fig. 5.3. We kunnen eveneens op een vast tijdstip (1.9 s) het temperatuursverloop in een snede beschouwen voor beide configuraties. We opteren voor de snede {x = 51, z = 51} wat een profiel oplevert in de y-richting langs de naalden. Dit wordt weergegeven in fig. 5.5.

zonder tumor met tumor

350

Temperatuur [K]

345 340 335 330 325 320 315 310 42

44

46

48

50

52

54

56

58

60

y [voxelnummers]

Figuur 5.5: Vergelijking van het temperatuursverloop in de snede {x = 51, z = 51} op het tijdstip 1.9 s in aan- of afwezigheid van een tumor.

Opnieuw merken we de verschillen tussen beide configuraties op. De temperatuur in de configuratie zonder tumor is groter dan die met tumor, wat eveneens naar voor komt in fig. 5.3. Bovendien kunnen we in fig. 5.5 opnieuw waarnemen dat de temperatuur vooral lokaal piekt ter hoogte van de naalden. Indien je je verder verwijdert dan 3 mm blijkt de temperatuur nauwelijks te verschillen van de basistemperatuur. Vervolgens kunnen we de invloed van de temperatuursafhankelijkheid van de elektrische geleidbaarheid (σ1 ) nagaan. We gebruiken de geometrie zoals beschreven in de vorige sectie, dus met goed en maligne weefsel. In de ene simulatie stellen we σ1 = 0 K−1 wat overeenstemt met een temperatuursonafhankelijke elektrische geleidbaarheid. In het andere geval gebruiken we de werkelijke σ1 = 0.015 K−1 . Indien het effect van σ1 = 0.015 K−1 klein of verwaarloosbaar zou zijn, dan moeten we voor de hele simulatieduur slechts één maal de potentiaal berekenen. Uit de Poisson-vergelijking volgt immers dat bij een temperatuursonafhankelijke σ de potentiaal φ onveranderd blijft. Indien het effect niet verwaarloosbaar is, moeten we voor elke tijdstap de potentiaal opnieuw berekenen zodat de rekentijd in dit geval veel groter zal zijn.

73


In fig. 5.6 wordt het temperatuursverloop ter hoogte van (51, 45, 51) tijdens één puls afgebeeld voor de twee verschillende waarden van σ1 . 355

σ1 = 0 K −1 σ1 = 0.015 K −1

350

Temperatuur [K]

345 340 335 330 325 320 315 310 0

1

2

3

4

5

6

7

Tijd [s]

Figuur 5.6: Vergelijking van het temperatuursverloop in het punt (51, 45, 51) voor σ1 = 0 K−1 en σ1 = 0.0015 K−1 .

Op basis van deze figuur moeten we concluderen dat σ1 wel degelijk een niet te verwaarlozen effect heeft op het temperatuursverloop. Helaas is het dus wel degelijk noodzakelijk om voor elke tijdstap een nieuw elektromagnetisch probleem te berekenen. Indien σ1 6= 0 K−1 wordt de maximumtemperatuur van 80 ◦ C sneller bereikt, zodat de stroom sneller wordt afgeschakeld en de temperatuur bijgevolg daalt. Ook spatiaal zijn er verschillen waar te nemen. Hiervoor bekijken we voor beide σ1 de snede z = 51 op die tijdstippen waarvoor elke configuratie zijn maximale temperatuur bereikt (tm,1 = 3.2 s voor σ1 = 0 K−1 en tm,2 = 2.7 s voor σ1 = 0.015 K−1 ). Het verschil in temperatuur tussen de beide gevallen wordt weergegeven in fig. 5.7 waarbij we voor de duidelijkheid enkel een gedeelte van de snede weergeven in de buurt van de bovenste naald. Deze figuur toont aan dat er een verschil van ongeveer 3.5 K bestaat tussen beide gevallen als ze hun maximum bereiken. Deze verschillen blijven echter beperkt tot op een afstand van ongeveer 4 mm. Dit valt gemakkelijk in te zien via het feit dat door het verschil in σ ook de potentialen zullen verschillen. Deze potentialen zijn enkel groot in de onmiddellijke nabijheid van de naald zodat het grootste verschil daar merkbaar zal zijn. Het sterk lokale karakter van de temperatuursverschillen komt eveneens naar boven in fig. 5.8. Daartoe bekijken we opnieuw de temperatuursverschillen tussen beide configuraties wanneer ze hun maximale temperatuur bereiken (tm,1 resp. tm,2 ). We berekenen dan de norm van de temperatuursverschillen binnen eenzelfde z-snede, de z-coördinaat ligt dus vast: zf . Deze norm wordt voor de verschillende zf afgebeeld in fig. 5.8. Uit deze figuur kunnen we zien dat, aangezien de naald stroom injecteert ter hoogte van zf = 51, de temperatuursverschillen slechts significant zijn tot zo’n 4 mm van de naaldtip.

74


68

3.5

66 3

y [voxelnummers]

64 2.5

62 60

2

58 1.5

56 54

1

52 0.5 50 48 35

40

45

50

55

60

65

x [voxelnummers]

Figuur 5.7: Spatiale weergave van het verschil in temperatuur op de tijdstippen van maximale temperatuur voor σ1 = 0 K−1 en σ1 = 0.0015 K−1 .

10

kT1 (x, y, zf , tm,1 ) − T2 (x, y, zf , tm,2 )k

9 8 7 6 5 4 3 2 1 0 46

48

50

52

54

56

58

zf -snede [voxelnummers]

Figuur 5.8: Norm van het temperatuursverschil binnen eenzelfde z-snede in functie van de co¨ ordinaat in de z-richting.

5.3

Invloed van discretisatie

Om een adequate ruimtelijke discretisatie te bepalen, gaan we de invloed van verschillende stapgroottes na. De referentiegeometrie bestaat uit 51×51×51 voxels met ∆x = ∆y = ∆z = 0.002 m. De straal van het goede weefsel (Rg ) bedraagt 20 voxels, terwijl de concentrische bol van maligne weefsel een straal heeft van 8 voxels (Rm ). De naalden worden ingebracht op (26, 30, 26) en (26, 22, 26) (roostercoördinaten). We bouwen nu identiek dezelfde omgevingen, maar met meer voxels bv. 75 × 75 × 75 of 101 × 101 × 101. Om eenzelfde geometrie te verkrijgen als deze met 51 × 51 × 51 voxels, moeten we dus ∆x = ∆y = ∆z en de andere geometrische objecten herschalen. Deze herschalingen worden samengevat in tabel 5.2.

75


Tabel 5.2: Afmetingen die de verschillende geometrieën karakteriseren om de invloed van discretisatie na te gaan.

Aantal voxels ∆x [m] Rg [voxels] Rm [voxels] Positie naald 1 [voxels] Positie naald 2 [voxels]

51 × 51 × 51 0.002 20 8 (26, 30, 26) (26, 22, 26)

75 × 75 × 75 51 0.002 · 75 30 12 (37, 43, 37) (37, 31, 37)

101 × 101 × 101 51 0.002 · 101 40 16 (51, 59, 51) (51, 44, 51)

We simuleren elk van deze geometrieën voor een tijdsspanne van twee pulsen. De benodigde rekentijd voor de uitvoering van de simulatie stijgt sterk met het aantal voxels zoals weergegeven in tabel 5.3. Tabel 5.3: Benodigde rekentijd voor de simulatie van de verschillende geometrieën.

Aantal voxels Rekentijd [min.]

51 × 51 × 51 5

75 × 75 × 75 12

101 × 101 × 101 65

Om de resultaten van de simulaties met elkaar in verband te brengen, bekijken we het temperatuursverloop langsheen de naalden in de y-richting op het tijdstip waar de maximale temperatuur behaald wordt. Dit wordt weergegeven in fig. 5.9. 370 ∆x = 0.002 ·

∆x = 0.002 m ∆x = 0.002 · 51 75 m

360

Temperatuur [K]

51 101 m

350

340

330

320

310 0.04

0.042

0.044

0.046

0.048

0.05

0.052

0.054

0.056

0.058

0.06

y [m]

Figuur 5.9: Temperatuursverloop in de y-richting langsheen de naalden op het tijdstip van maximale temperatuur voor de verschillende discretisatieniveaus.

We merken op dat de vorm van de pieken in de drie gevallen vrij gelijklopend is. De pieken zijn een beetje verschoven ten opzichte van elkaar door discretisatiefouten die ervoor zorgen dat de simulatieomgeving niet exact gelijk is. Bij het geval voor ∆x = 0.002 m is het duidelijk dat de stapgrootte te ruw is zodat de maximale temperatuur significant hoger ligt dan 80 ◦ C. De

76


afmetingen van dergelijke voxels zijn immers vrij groot zodat de modellering van de diffusie en het elektromagnetisch deelprobleem minder nauwkeurig uitgevoerd wordt en de temperatuur 51 gedurende één tijdstap (0.1 s) snel kan stijgen. Het is pas bij ∆x = 0.002 · 101 m ≈ 0.001 m dat het beslissingscriterium nauwkeurig toegepast kan worden. Wij opteren altijd voor een discretisatie met ∆x = 0.001 m omdat we met behulp van een dergelijk discretisatieniveau in een aanvaardbare rekentijd mooie simulaties kunnen creëren. We moeten dus niet naar resoluties gaan op µm-schaal.

5.4

Interpretatie van de numerieke simulatie

We kunnen de numerieke simulaties ook aangrijpen om de invloed van andere effecten te interpreteren. Zo kunnen we eens kijken wat de invloed is van de bloedperfusie in het haarvatennetwerk van de lever. Aangezien de massadichtheid en de specifieke warmtecapaciteit constant zijn, komt dit overeen met een studie van ω. Deze constante stelt de hoeveelheid bloed voor die per volume-eenheid en per seconde door het weefsel stroomt. Uit tabel 3.2 volgt dat ω = 0.0064 s−1 . Indien ω = 0 s−1 zou zijn, dan komt dit in de realiteit overeen met de situatie waarin er geen bloeddoorstroming is. Dit kan tijdens een RFA-procedure bekomen worden door de toevoerslagaders af te klemmen. We simuleren de beide gevallen voor ω in de simulatieomgeving beschreven in sectie 5.1 gedurende één cyclus. In fig. 5.10 bekijken we het temperatuursverloop in het punt (51, 45, 51) in de onmiddellijke nabijheid van een naaldtip. 355

ω = 0 s−1 ω = 0.0064 s−1

350

Temperatuur [K]

345 340 335 330 325 320 315 310 0

1

2

3

4

5

6

Tijd [s]

Figuur 5.10: Vergelijking van het temperatuursverloop in het punt (51, 45, 51) met of zonder bloedperfusie.

Op basis van deze figuren zouden we kunnen concluderen dat de bloedperfusie slechts een marginaal effect uitoefent op het temperatuursverloop. Het hoeft echter niet te verwonderen dat de bloedperfusie weinig 2 effect heeft ter hoogte van de naalden. Op die posities zal de ~ elektrische bronterm σ E immers dominant zijn. Op plaatsen die verder verwijderd zijn van

de naald zal deze bronterm veel minder groot zijn, zodat de bloedperfusie wel een belangrijk effect zal hebben. Zo kunnen we de norm beschouwen van het verschil in temperatuur in de

77


snede z = 53 tussen de twee gevallen voor ω. Deze norm wordt weergegeven in fig. 5.11 als functie van de tijd. 0.5

zf =53

kT1 (x, y, zf , t) − T2 (x, y, zf , t)k

0.45 0.4

0.35 0.3

0.25 0.2

0.15 0.1

0.05 0 0

1

2

3

4

5

6

t [s]

Figuur 5.11: Norm van de temperatuursverschillen in de snede zf = 53 in functie van de tijd.

Hier verkrijgen we wel degelijk een observeerbaar verschil in temperatuur. Bovendien stijgt de norm van de temperatuursverschillen monotoon in de tijd. Hoe langer we de RFA-procedure uitvoeren, hoe groter het verschil tussen de twee simulaties. Merk op dat dit enkel en alleen door het haarvatennetwerk in de lever wordt veroorzaakt. In de lever bevinden zich ook tal van grotere bloedvaten die ongetwijfeld een groot effect zullen hebben op de temperatuursverdeling in de lever. In subsectie 2.3.1.3 en in de literatuur werd dit heat sink effect vermeld. Dit effect kunnen we met ons model niet inschatten. Grotere bloedvaten zijn immers niet zo eenvoudig te modelleren in een continu¨ ummodel. Een eventuele oplossing zou kunnen zijn om een convectieterm met specifieke parameters in te bouwen in de Pennes bio-warmtevergelijking. Dit is echter buiten het bereik van deze masterproef. Verder kunnen we ook de invloed analyseren van onze modulatie van de Pennes bio-warmtevergelijking zoals voorgesteld in sectie 3.6. De verwachting is dat deze modulatie vooral een effect zal hebben bij het afschakelen van de stroom. Dan is de bloedperfusieterm een belangrijke term. Aangezien 0 ≤ α ≤ 1 zal de warmteafvoerende capaciteit van de haarvaten verminderd worden en de tijdsconstante van het temperatuursverval groter zijn. Bovendien zal het verschil tussen de simulatie met of zonder modulatie alleen maar toenemen gezien het feit dat de overlevingsgraad monotoon dalend is. We simuleren de situatie met en zonder modulatie gedurende twee pulsen. Het temperatuursverloop in (51, 45, 51) wordt weergegeven in fig. 5.12. In deze figuur is duidelijk zichtbaar dat er wel degelijk verschillen zijn tussen beide simulaties. Gedurende de eerste aanschakelperiode lopen beide oplossingen gelijk aangezien het weefsel nog weinig schade heeft geleden en dus α ≈ 1. Desalniettemin zijn er kleine verschillen zodat de configuratie met de modulatie iets vroeger de maximale temperatuur bereikt en dus een temperatuursval ondergaat. Gedurende de afschakeltijd wordt het verschil tussen de beide temperaturen steeds kleiner zodat we mogen oordelen dat de configuratie met de modula-

78


355

zonder modulatie met modulatie

350

Temperatuur [K]

345 340 335 330 325 320 315 310 0

2

4

6

8

10

12

Tijd [s]

Figuur 5.12: Vergelijking van het temperatuursverloop in het punt (51, 45, 51) met of zonder modulatie van de Pennes bio-warmtevergelijking.

tie inderdaad trager vervalt. Tijdens de tweede maal dat de stroom aangeschakeld wordt, zien we het tegenovergestelde. De simulatie met de modulatie stijgt sneller dan degene zonder modulatie. Opnieuw is dit te verwachten door het feit dat de simulatie met modulatie rekening houdt met de gediminueerde warmteafvoerende capaciteit. Gedurende de tweede afschakeltijd neemt het verschil in temperatuur weer af zodat de oplossing met de modulatie trager vervalt. We kunnen ook onmiddellijk kijken naar de overlevingsgraad zelf (fig. 5.13). 1

zonder modulatie met modulatie

Overlevingsgraad [-]

0.9 0.8 0.7 0.6 0.5 0.4

0

2

4

6

8

10

12

Tijd [s]

Figuur 5.13: Vergelijking van de overlevingsgraad in het punt (51, 45, 51) met of zonder modulatie van de Pennes bio-warmtevergelijking.

Als de stroom aanstaat, daalt de overlevingsgraad relatief sterk aangezien er schade wordt aangericht. Door het feit dat de stroom in het geval met modulatie iets sneller wordt stopgezet (zie fig. 5.12), zal zijn overlevingsgraad rond dat tijdstip ook groter zijn. Er is immers iets


79

minder schade aangericht. Tijdens de situatie wanneer de stroom afgeschakeld is, merken we slechts kleine veranderingen in overlevingsgraad op. Dit is logisch aangezien de temperatuur in dit gebied vrij laag ligt, zodat er weinig schade wordt aangericht. Pas vanaf het moment dat de temperatuur weer significant wordt, daalt de overlevingsgraad opnieuw. Tijdens de tweede aanschakeltijd is opnieuw merkbaar dat de modulatie de temperatuur sterker doet toenemen. De overlevingsgraad bij het begin van het aanschakelen is immers groter in het geval van de modulatie, terwijl de overlevingsgraad bij het moment van afschakelen (7.6 s) in de beide gevallen dezelfde is. Er is dus tijdens het aanschakelen meer schade aangericht in de situatie waar we de modulatie toepassen. Ten slotte kunnen we uit de simulaties in dit hoofdstuk twee belangrijke argumenten halen die het succes van RFA ernstig beperken. Vooreerst hebben we in verschillende onafhankelijke simulaties kunnen waarnemen dat de elektrische potentialen, temperatuursverhogingen en overlevingsgraad slechts veranderen in de onmiddellijke nabijheid van de naaldtippen. Afhankelijk van simulatie tot simulatie beperkt de invloedssfeer van RFA zich tot een regio met straal 3 ` a 4 mm. Uiteraard simuleren we in dit hoofdstuk slechts een heel korte RFA-procedure. Naarmate er langer stroom geleverd wordt, zal de conductie van warmte ervoor zorgen dat de invloedssfeer zal uitbreiden. Desalniettemin moeten we concluderen dat de toepassing van RFA resulteert in een heel lokale behandeling. Dit is enerzijds voordelig aangezien er geen schade kan aangericht worden in grote delen van het gezonde weefsel. Anderzijds wordt het zo moeilijk om zonder al te veel inserties grotere tumoren te behandelen. Een tweede beperking ligt in de invloed van de grote bloedvaten. Wij bestudeerden de bloedperfusie van het haarvatennetwerk dat een niet te verwaarlozen effect heeft op de temperatuursdistributie. Grotere bloedvaten veroorzaken via convectie een heat sink effect. Dit belemmert het uitvoeren van een succesvolle RFA. In deze masterproef spenderen we geen aandacht aan het tweede effect. In hoofdstuk 6 zullen we trachten het eerste effect te beperken via een optimalisatie van de stroomvorm.

Hoofdstuk 6

Verbeterde thermotherapie voor leverkanker There’s a way to do it better - find it. Thomas Edison

6.1

Inleiding

Zoals aangehaald op het einde van het vorige hoofdstuk zijn er beperkingen bij het toepassen van thermotherapie onder de vorm van RFA. Een van de grootste voordelen, het lokaal zijn van de thermotherapie, blijkt ook een fundamenteel nadeel te zijn. Elke naald kan slechts een object met een beperkte straal aanvallen. Zo blijkt thermotherapie niet in staat te zijn om grotere tumoren succesvol te behandelen. Tegenwoordig gebruikt men de monopolaire configuratie voor RFA. Hierbij gebruikt men één enkele naald die één sferische ablatiezone kan veroorzaken. Ook in deze configuratie zijn er grenzen aan de grootte van de tumoren. Om grotere tumoren te behandelen moet men meerdere inserties gebruiken tijdens de RFA-procedure of moet men zijn toevlucht nemen tot multitine-elektrodes waarbij de naald eindigt op verschillende tippen in een parapluvorm. Zo kan men grotere volumes trachten te behandelen. In deze masterproef bestuderen we de bipolaire configuratie. Hier werken we met naalden die door een stroomkring aangestuurd worden. Zo verkrijgen we twee ablatiezones. Elk van deze zones is in staat om een klein volume tumorweefsel te omvatten. Het grote voordeel is echter dat die zones tegelijkertijd ontstaan en er dus een mogelijkheid bestaat dat de ene ablatiezone de andere be¨ınvloedt. Door de overlap kunnen we een versterking van de temperatuur bekomen op het scheidingsoppervlak van de ablatiezones wat kan resulteren in een efficiëntere behandeling van grotere tumoren. Deze overlap zal in belangrijke mate be¨ınvloed worden door het ge¨ınjecteerde stroomprofiel. In dit hoofdstuk zullen we dit concept uitwerken.

80

81

Hoofdstuk 6. Verbeterde thermotherapie voor leverkanker

6.2 6.2.1

Gepulste thermotherapie Simulatieomgeving

In dit hoofdstuk zullen we een omgeving beschouwen bestaande uit 51 × 51 × 51 voxels met ∆x = ∆y = ∆z = 0.001 m. In deze omgeving plaatsen we een sferisch volume gezond weefsel met een straal Rg = 20 voxels. Binnen het gezonde weefsel plaatsen we een concentrisch volume maligne weefsel met een straal Rm = 8 voxels. De naaldtippen worden ingebracht ter hoogte van (26, 30, 26) en (26, 22, 26) zodat de afstand van elke naald tot het centrum (26, 26, 26) 4 mm bedraagt. In deze configuratie zal de invloed van een verschillende stroomvorm het duidelijkst naar voor komen.

6.2.2

Referentiestroomvorm

In de literatuur besteedt men weinig aandacht aan de stroomvorm bij een RFA-procedure. We zijn dus genoodzaakt om een eigen referentiestroomvorm te bepalen. In sectie 5.1 motiveerden we ons initieel beslissingscriterium inzake de temperatuur. Onze referentiestroomvorm wordt grafisch weergegeven in fig. 6.1. I / I 0 - 1

Tmax t1   80 C

Tmax t2   50 C

t1

t2

t s 

Figuur 6.1: Referentiestroomvorm corresponderend met het beslissingscriterium inzake de temperatuur.

Om te vermijden dat er onvoorziene effecten ontstaan, schakelen we de stroom af wanneer de maximumtemperatuur in het weefsel 80 ◦ C bedraagt. Vervolgens laten we de temperatuur relaxeren tot de maximale temperatuur gezakt is tot 50 ◦ C waar er quasi geen thermische schade meer wordt aangericht. Door het feit dat er een boven- en een benedengrens is voor de maximumtemperatuur resulteert het stroomprofiel in een gepulseerde injectie van stroom. 6.2.2.1

Numerieke simulatie

We simuleren een RFA-procedure van 10 pulsen met een referentiestroomvorm. Figuur 6.2 geeft een visuele weergave van het ge¨ınjecteerde stroomprofiel. Op deze figuur kunnen we de tien pulsen waarnemen. Bij elk tijdstip waarop de stroom wordt

82


0.07

1

4

2

6

8

12

10

16

14

20

18

0.06

Stroom [A]

0.05 0.04 0.03 0.02 0.01 3 5

0 0

7 15

5 10

9 20

13 30

11 25

15 35

21

19

17 40

45

50

Tijd [s]

Figuur 6.2: Stroomprofiel gedurende een RFA-procedure van 10 pulsen met aanduiding van de pulslabels.

aan- of afgeschakeld vermelden we een getal, het pulslabel. Dit pulslabel zal gebruikt worden als ordinaat om de temperatuur of overlevingsgraad uit te zetten. Als resultaat van de simulatie geven we in fig. 6.3 de overlevingsgraad op verschillende locaties weer in functie van de pulslabels. 1

nabij centrum

0.9


0.8 0.7

dichtbij naald

0.6 0.5 0.4 0.3 0.2 0.1 0 1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

Pulslabel [-]

Figuur 6.3: Overlevingsgraad op locaties dichtbij de naald en het centrum in functie van de pulslabels.

Dichtbij de naald daalt de overlevingsgraad vrij sterk zodat het omringende maligne weefsel na 10 pulsen al vernietigd is. Dit toont aan dat RFA heel doeltreffend is op korte afstand. Dichtbij het centrum op de grens van de twee ablatiezones zien we echter een ander beeld. De overlevingsgraad daalt hier traag zodat we na 10 pulsen maar weinig schade hebben aangericht. In het vervolg van dit hoofdstuk zullen we met bijzondere interesse uitkijken naar de temperatuur en overlevingsgraad ter hoogte van het centrum (26, 26, 26). Dit centrum

83


is immers gesitueerd in het midden van de verbindingslijn tussen de naaldtippen. In die hoedanigheid is dit punt uitermate geschikt om de overlap te karakteriseren tussen beide ablatiezones. In fig. 6.4 geven we de temperatuur na elke puls weer ter hoogte van het centrum. 318 317

Temperatuur [K]

316 315 314 313 312 311 310 0

1

2

3

4

5

6

7

8

9

10

Aantal pulsen [-]

Figuur 6.4: Temperatuur na elke puls ter hoogte van het centrum (26, 26, 26).

Het is duidelijk dat de temperatuur in het centrum toeneemt met elke puls. Het dominerende mechanisme dat zorgt voor deze temperatuursverhoging is de conductie van warmte uit de zones in de buurt van de naalden. Hoe meer pulsen we aanleggen, hoe hoger de temperatuur zal zijn in dit centrum en hoe meer schade we kunnen aanrichten.

6.3 6.3.1

Verbeterde thermotherapie Concept voor een effici¨ entere thermotherapie

De aanleiding om de variatie van de stroomvorm voorop te stellen in de zoektocht naar een verbeterde thermotherapie wordt teruggevonden in fig. 6.5. Deze figuur correspondeert met vgl. (4.73). De startsituatie betreft een beperkte zone met een temperatuur van 100 K, terwijl de temperatuur in de overige gebieden 0 K bedraagt. Vervolgens geeft vgl. (4.73) het temperatuursverval weer in de afwezigheid van stroom. Fig. 6.5 geeft op verschillende tijdstippen het temperatuursprofiel weer in de centrale snede {y = 26, z = 26}. Op basis van deze figuur kunnen we twee zaken opmerken. Enerzijds zakt de maximumtemperatuur naarmate de tijd vordert, anderzijds zal het temperatuursprofiel breder worden. Deze observatie kunnen we linken met onze zoektocht naar een verbeterde thermotherapie. Om grotere tumoren te kunnen behandelen willen we net dat er schade wordt aangericht in een grotere zone rond de naaldtippen. Daartoe moeten we vooral de invloed van de afschakeltijd τ analyseren. Het betreft een typisch ingenieursprobleem met tegengestelde invloeden. Als we de stroom maar een korte periode afschakelen, zal de maximumtemperatuur groot zijn.

84


100 90

Temperatuur [K]

80 70 60 50 40 30 20 10 0 −0.025

−0.02

−0.015

−0.01

−0.005

0

0.005

0.01

0.015

0.02

0.025

x [m]

Figuur 6.5: Temperatuursverval in de tijd indien geen stroom geleverd wordt.

Helaas zal de breedte van het temperatuursprofiel nagenoeg niet toegenomen zijn. We zullen dus grote schade aanrichten in een klein volume, terwijl de overige regio’s quasi onaangeroerd blijven. Schakelen we de stroom gedurende een lange tijdsspanne af, verkrijgen we een heel breed temperatuursprofiel. De maximumtemperatuur zal echter zo laag liggen dat er geen thermische schade zal aangericht worden over het hele profiel. We hopen een significante verbetering waar te nemen bij een intermediaire afschakeltijd. We kunnen de afschakeltijd via twee verschillende beslissingscriteria be¨ınvloeden. Enerzijds kunnen we opteren om de stroom opnieuw aan te schakelen bij een hogere temperatuur dan 50 ◦ C zoals bij onze referentiestroomvorm. Anderzijds kunnen we een vaste afschakeltijd τf opleggen. De stroom zal dan gedurende een vaste tijdsperiode afgeschakeld zijn onafhankelijk van de maximale temperatuur. We zullen beide stroomvormen in de volgende subsecties analyseren.

6.3.2

Aanschakeltemperatuur

Vooreerst bekijken we het effect van een andere aanschakeltemperatuur. We injecteren een analoge stroomvorm als de referentiestroomvorm uit fig. 6.1, alleen veranderen we de temperatuur van 50 ◦ C om de stroom weer aan te schakelen. In fig. 6.6 geven de temperatuur weer in het centrum (26, 26, 26) voor verschillende aanschakeltemperaturen. De temperatuur wordt hier uitgezet tegenover het aantal pulsen dat we aangelegd hebben. Zoals duidelijk weergegeven blijkt de temperatuur sterk af te hangen van de aanschakeltemperatuur. Na circa 5 pulsen geeft de referentiestroomvorm met T = 50 ◦ C de hoogste temperatuur in het centrum. In het gebied van 20 tot 50 pulsen zorgen T = 65 ◦ C en T = 70 ◦ C echter voor de hoogste temperaturen. Ook in de overlevingsgraad vinden we een sterke afhankelijkheid terug van de aanschakeltemperatuur (fig. 6.7). Tot pulslabel 30 blijkt de referentiestroomvorm de grootste schade aan te richten. Dit is uiteraard in overeenstemming met het feit dat deze stroomvorm de hoogste temperatuur be-

85


324 322

Temperatuur [K]

320 318 316

T T T T T

314 312 310 0

5

10

15

20

25

30

35

40

= 50 ◦ C = 55 ◦ C = 60 ◦ C = 65 ◦ C = 70 ◦ C

45

50

Aantal pulsen [-]

Figuur 6.6: Temperatuur na elke puls ter hoogte van het centrum voor verschillende aanschakeltemperaturen.

1

T T T T T

0.9

= 50 ◦◦ C = 55 ◦ C = 60 ◦ C = 65 ◦ C = 70 C


0.8 0.7 0.6 0.5 0.4 0.3 0.2 0

20

40

60

80

100

120

Pulslabel [-]

Figuur 6.7: Overlevingsgraad ter hoogte van het centrum bij elk pulslabel en voor verschillende aanschakeltemperaturen.

reikt in het begin zoals weergegeven in fig. 6.6. Vanaf pulslabel 30 tot pulslabel 80 levert de stroomvorm met aanschakeltemperatuur T = 55 ◦ C de beste resultaten op. Na 50 pulsen (pulslabel 101) blijkt echter dat een aanschakeltemperatuur van T = 60 ◦ C de meeste schade heeft aangericht. Zoals vermeld bij de vorige figuur bereikt de aanschakeltemperatuur T = 70 ◦ C de grootste temperatuur in het gebied van 20 tot 50 pulsen. Dit is eveneens observeerbaar in fig. 6.7. In dit gebied zakt de overlevingsgraad van de simulatie met deze

86


aanschakeltemperatuur immers het snelst. Bovendien is het zo dat een simulatie met aanschakeltemperatuur T = 70 ◦ C minder stroom zal leveren gedurende een puls ten opzichte van de referentiestroomvorm. Hoe hoger de aanschakeltemperatuur, hoe minder tijd er immers nodig is om de maximumtemperatuur van T = 80 ◦ C te bereiken. Het aantal seconden dat er stroom wordt geleverd verschilt sterk tussen de verschillende aanschakeltemperaturen. Zo is voor een aanschakeltemperatuur van T = 70 ◦ C slechts gedurende 36.8 s stroom geleverd na 50 pulsen, terwijl de referentiestroomvorm gedurende 81.6 s stroom heeft ge¨ınjecteerd. Willen we analyseren welke aanschakeltemperatuur de meest efficiënte therapie oplevert, dan moeten we hiermee rekening houden. In fig. 6.8 geven we de overlevingsgraad weer in functie van het aantal seconden dat er stroom geleverd wordt. 1

T T T T T

0.9

= 50 ◦◦ C = 55 ◦ C = 60 ◦ C = 65 ◦ C = 70 C


0.8 0.7 0.6 0.5 0.4 0.3 0.2 0

10

20

30

40

50

60

70

80

90

Aantal seconden stroom [s]

Figuur 6.8: Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende aanschakeltemperaturen.

Verrassend genoeg moeten we vaststellen dat het duidelijk is dat de stroomvorm met aanschakeltemperatuur T = 70 ◦ C veruit het efficiëntst is. Hoe hoger de aanschakeltemperatuur, hoe meer schade er wordt aangericht bij eenzelfde aantal seconden dat er stroom wordt ge¨ınjecteerd. Dit resultaat is uiteraard afhankelijk de afstand van de naalden tot het centrum. Bij onze simulatieomgeving is er slechts een afstand van 4 mm tussen de beide punten, zodat het diffusieproces van warmte relatief snel kan plaatsgrijpen.

6.3.3

Vaste afschakeltijd

Een tweede mogelijkheid om de stroomvorm aan te passen is door niet langer de afschakeltijd te laten afhangen van de maximumtemperatuur. In plaats van de stroom opnieuw aan te schakelen als de maximumtemperatuur gezakt is tot bv. 50 ◦ C, leggen we een vaste afschakeltijd τf op. Gedurende de tijdsspanne τf injecteren we geen stroom meer, waarna de cyclus opnieuw begint. We bekijken de overlevingsgraad in het centrum voor verschillende τf in fig. 6.9.

87


1


0.9


0.8 0.7

= 0.5 s = 1.0 s = 1.5 s = 2.0 s = 2.5 s = 3.0 s

0.6 0.5 0.4 0.3 0.2 0

20

40

60

80

100

120

Pulslabel [-]

Figuur 6.9: Overlevingsgraad ter hoogte van het centrum bij elk pulslabel en voor verschillende afschakeltijden.

Deze figuur toont opnieuw aan dat de afschakeltijd een grote invloed heeft op de overlevingsgraad. Bij elk pulslabel geeft een korte afschakeltijd van 0.5 s op het eerste zicht de slechtste resultaten. Over de volledige simulatieduur blijken intermediaire afschakeltijden (τf = 1.5, 2.0, 2.5 s) de beste resultaten op te leveren. Vanaf 30 pulsen (pulslabel 61) wordt de meeste schade aangericht door de afschakeltijd τf = 2.0 s. Voor deze afschakeltijd geven we in fig. 6.10 eveneens een spatiale weergave van de overlevingsgraad weer. 50

start

40

y [voxelnummers]

y [voxelnummers]

50

30 20

10de puls 1

40

0.9

30

0.8

20

0.7 10

10

0.6 10 50

20

30

x [voxelnummers]

40

10

50 50

20ste puls

20

30

x [voxelnummers]

40

50 0.5

50ste puls

y [voxelnummers]

y [voxelnummers]

0.4 40 30 20 10

10

20

30

x [voxelnummers]

40

50

40 0.3 30

0.2

20

0.1

10

0

10

20

30

x [voxelnummers]

40

50

Figuur 6.10: Spatiale weergave van de overlevingsgraad in de snede z = 26.

88


Na tien pulsen zijn er twee afzonderlijk ablatiezones aanwezig met een straal van circa 2.5 mm. Ter hoogte van het centrum is de temperatuur te laag waardoor er nog geen schade is aangericht. Naarmate er meer pulsen aangelegd worden, zal de temperatuur in dit centrum stijgen en de overlevingsgraad dalen. Na 50 pulsen is het centrum wel degelijk aangetast door de overlap van de twee ablatiezones. Door enkel en alleen te spelen met de afschakeltijd kunnen we grotere objecten behandelen. Om de vraag te kunnen beantwoorden welke afschakeltijd nu de efficiëntste therapie oplevert, moeten we zoals in de vorige subsectie opnieuw rekening houden met de tijd waarin er stroom geleverd wordt. Immers, na een korte afschakeltijd zal de maximumtemperatuur nog niet zo sterk afgenomen zijn. In een volgende puls zal er dus minder tijd nodig zijn om een maximumtemperatuur van 80 ◦ C. Zo wordt er voor 50 pulsen slechts gedurende 29.7 s stroom ge¨ınjecteerd bij τf = 0.5 s tegenover 73.5 s voor τf = 3.0 s. Figuur 6.11 geeft een visuele weergave van de temperatuur in het centrum in functie van het aantal seconden dat er stroom ge¨ınjecteerd wordt. 324 322

Temperatuur [K]

320 318 316


314 312 310 0

10

20

30

40

50

60

= 0.5 s = 1.0 s = 1.5 s = 2.0 s = 2.5 s = 3.0 s

70

80


Figuur 6.11: Temperatuur ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden.

Indien we rekening houden met het aantal seconden dat er stroom gevoerd wordt, dan zien we dat er geen verschil in temperatuur is gedurende de eerste 10 s. Vanaf dan bereiken we een hogere temperatuur in het centrum indien we een afschakeltijd van 0.5 s gebruiken. Hoe korter τf , des te hoger de temperatuur in het centrum bij eenzelfde tijdsduur van stroominjectie. Bij de overlevingsgraad merken we hetzelfde fenomeen op (fig. 6.12). Door het feit dat de temperatuur hoger is voor τf = 0.5 s, is ook de overlevingsgraad veel lager in dit geval. Deze korte afschakeltijd blijkt dus de meest efficiënte vorm van RFA voort te brengen in deze simulatieomgeving. Dit is op zijn minst merkwaardig te noemen in het licht van fig. 6.9. In functie van het aantal pulsen bleek τf = 0.5 s nog de slechtste resultaten op te leveren. Het feit dat de tijdsduur van elke puls zo kort is voor τf = 0.5 s zorgt voor dit vertekend beeld.

89


1


0.9


0.8

= 0.5 s = 1.0 s = 1.5 s = 2.0 s = 2.5 s = 3.0 s

0.7 0.6 0.5 0.4 0.3 0.2 0

10

20

30

40

50

60

70

80


Figuur 6.12: Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden.

We kunnen met de optimale afschakeltijd voor deze configuratie, nl. τf = 0.5 s, een langere simulatie uitvoeren. We opteren voor een simulatie met 500 pulsen zodat er met deze afschakeltijd gedurende 228.1 s stroom wordt ge¨ınjecteerd. Vooreerst bekijken we in fig. 6.13 opnieuw de overlevingsgraad in het centrum als maatstaf voor de overlap van ablatiezones. 1 0.9


0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

50

100

150

200


Figuur 6.13: Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom.

Uiteraard zal er in de eerste 10 s weinig schade aangericht worden aangezien de temperatuur in de centrale voxel nog niet verhoogd is. Vanaf het tijdstip dat de temperatuur toeneemt,

90


neemt de overlevingsgraad sterk af. Na slechts 60 s stroom is het tumorweefsel in het centrum volledig vernietigd. Uiteraard heeft deze korte tijd deels te maken met de kleine afmetingen van de simulatieomgeving, desalniettemin kunnen we ook de invloed van de optimale afschakeltijd aanhalen. Zoals zichtbaar in fig. (6.12) zijn de verschillen tussen de verschillende afschakeltijden al vrij groot na 50 pulsen. Na 500 pulsen wordt dit verschil alleen maar groter zodat een optimalisatie voor elke configuratie wel degelijk zinvol is. Om het uiteindelijke resultaat van deze simulatie te analyseren, bekijken we in fig. 6.14, 6.15 en 6.16 spatiale weergaves van de overlevingsgraad in de centrale sneden op het eindtijdstip.

50 0.9

45

0.8

z [voxelnummers]

40

0.7

35

0.6

30

0.5

25

0.4

20 15

0.3

10

0.2

5

0.1 10

20

30

40

50

0

y [voxelnummers]

Figuur 6.14: Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede x = 26.

In fig. 6.14 (de snede x = 26) vinden we een rechthoekig ablatiegebied terug met basis 7 mm en een hoogte die 12 mm bedraagt. Op basis van deze figuur kunnen we bepalen dat de afzonderlijke ablatiezones van de twee naalden volledig samengesmolten zijn tot één grote zone. Ter hoogte van het centrum vindt er dus een temperatuursopbouw plaats die voldoende hoog is om de tumorcellen te doden. Een interessant beeld vormt de visuele weergave in fig. 6.15. Hier wordt de overlevingsgraad weergegeven in het middelloodvlak van de naalden die in de y-richting ingebracht worden. De ablatiezone in de snede heeft de vorm van een cirkel met diameter 7 mm. De circulaire vorm van de individuele ablatiezones blijft dus behouden in de globale ablatiezone. Kijken we ten slotte in fig. 6.16 naar de overlevingsgraad in de snede z = 26, dan zien we opnieuw het beeld van fig. 6.14: een rechthoekig ablatiegebied. Uit deze drie spatiale weergaves van de overlevingsgraad kunnen we concluderen dat de bipolaire configuratie bij RFA een cilindrische ablatiezone veroorzaakt. De individuele sferische ablatiezones ter hoogte van elke naald kunnen met behulp van een optimale afschakeltijd samensmelten en het tussenliggende gebied volledig vernietigen. Zo kunnen we nog steeds gebruik maken van de lokale eigenschappen van RFA, maar blijven de ablatiezones niet beperkt tot een aantal mm in de buurt van de naalden. Dit principe zou kunnen helpen in de

91


50 0.9

45

0.8

z [voxelnummers]

40

0.7

35

0.6

30

0.5

25

0.4

20 15

0.3

10

0.2

5

0.1 10

20

30

40

50

0

x [voxelnummers]

Figuur 6.15: Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede y = 26. 50 0.9

45

0.8

y [voxelnummers]

40

0.7

35

0.6

30

0.5

25

0.4

20 15

0.3

10

0.2

5

0.1 10

20

30

40

50

0

x [voxelnummers]

Figuur 6.16: Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede z = 26.

toepassing van RFA in de behandeling van grotere tumoren. We wensen nogmaals uitdrukkelijk te vermelden dat de vaststelling dat τf = 0.5 s de meest efficiënte afschakeltijd is, afhankelijk is van de locatie van de naalden. Indien we de naalden verder van elkaar verwijderen, zal een intermediaire afschakeltijd het efficiëntst zijn. In fig. 6.17 wordt de overlevingsgraad weergegeven als we de naalden op de posities (26, 33, 26) en (26, 19, 26) inbrengen. Hierdoor bedraagt de afstand van elke naald tot het centrum 7 mm. Uit deze figuur kunnen we afleiden dat τf = 3.0 s de meest efficiënte afschakeltijd is. Ondanks het feit dat we rekening houden met het aantal seconden dat er stroom wordt geleverd, is de kortste afschakeltijd dus niet altijd de beste. Verder kunnen we opmerken dat hoe verder de naalden verwijderd zijn van het centrum, hoe minder schade er aangericht is in het cen-

92


1.002


1


0.998

= 0.5 s = 1.0 s = 1.5 s = 2.0 s = 2.5 s = 3.0 s

0.996 0.994 0.992 0.99 0.988 0.986 0.984 0

10

20

30

40

50

60

70


Figuur 6.17: Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden.

trum. In fig. 6.17 bedraagt de overlevingsgraad steeds meer dan 0.986 wat impliceert dat er quasi geen schade is aangericht. Om toch schade aan te richten, moeten we een langere RFA-procedure initialiseren met meer pulsen. Uit dit alles kunnen we concluderen dat de stroomvorm een grote invloed heeft in de behandeling van grotere tumoren. Bovendien laten de numerieke methoden ontwikkeld in het kader van deze masterproef toe om voor elke patiënt of tumorgrootte na te gaan wat de optimale stroomvorm is. Deze stroomvorm zal afhankelijk zijn van de afschakeltijd en de locatie van de naalden. In principe moeten we dus de therapie optimaliseren door voor elk mogelijk koppel van naaldposities het effect van de verschillende afschakeltijden te analyseren. Om te bepalen welke therapieopstelling het beste is, kunnen we de statistische verwerking gebruiken zoals uiteengezet in subsectie 3.5.2. Vooraleer we overgaan tot een dergelijke rekenintensieve analyse, is er eerst een experimentele validatie nodig. In deze masterproef beperken we ons echter tot in silico-experimenten.

Hoofdstuk 7

Conclusie I think and think for months and years, ninety-nine times, the conclusion is false. The hundredth time I am right. Albert Einstein

In deze masterproef hebben we getracht om de volgende centrale vragen te beantwoorden: “Hoe kan een stabiele numerieke methode ontwikkeld worden voor de thermotherapie van leverkanker?”, “Wat is de invloed van materiaalparameters op de spatio-temporele temperatuursdistributie?”, “Hoe kan de huidige techniek verbeterd worden voor een meer efficiënte behandeling van levertumoren?”

7.1

Overzicht van de masterproef

In hoofdstuk 1 formuleren we de probleemstelling en de daaraan gekoppelde doelstellingen van deze masterproef. Het werkingsprincipe van RFA wordt door middel van een literatuurstudie uiteengezet in hoofdstuk 2. Belangrijke elementen in deze literatuurstudie betreffen de wiskundige modellering van de biologische warmteoverdracht en de materiaaleigenschappen van gezonde en maligne levercellen. Door de biologische warmteoverdracht te bestuderen, bekwamen we de Pennes bio-warmtevergelijking. Deze vergelijking geeft aan hoe de temperatuur in een biologisch weefsel verandert in ruimte en tijd als gevolg van warmtebronnen en -afvoeren. Dit vormt de basis van het thermische probleem. Aan het einde van dit hoofdstuk bespreken we eveneens de monopolaire en bipolaire configuratie voor RFA. In deze masterproef bestuderen we de bipolaire configuratie omdat we in deze configuratie meer vrijheid hebben om te zoeken naar verbeteringen. Hoofdstuk 3 omvat de wiskundige modellering en de numerieke methoden voor het elektromagnetisch en thermisch deelprobleem. Verder schenken we eveneens aandacht aan de interpretatie van de resultaten wat resulteert in de overlevingsgraad op basis van een Arrheniusmodel voor de thermische schade. Op basis van deze overlevingsgraad stellen we een modulatie van de huidige Pennes bio-warmtevergelijking voor. Eveneens leggen we in dit hoofdstuk de theoretische basis om tot een stabiele numerieke methode te komen. Zo worden er voorwaarden geformuleerd op de temporele staplengten met behulp van de Fourier methode van Von Neumann. Hierna voeren we in hoofdstuk 4 een numerieke analyse door van onze 93

Hoofdstuk 7. Conclusie

94

gehele numerieke methode. De klemtoon ligt hierbij op het valideren van de verschillende deelproblemen. Analytische limietsituaties worden daartoe berekend om een vergelijking met de numerieke resultaten toe te laten. Eveneens gaan we de nauwkeurigheid na van het stabiliteitscriterium afgeleid in hoofdstuk 3. Op basis van de resultaten in dit hoofdstuk mogen we concluderen dat we voldoen aan het eerste doel van de thesis: de ontwikkeling van een stabiele numerieke methode. De simulatie van thermotherapie wordt beschreven in hoofdstuk 5. Hierbij hebben we een antwoord geformuleerd op de vraag naar de invloed van de materiaalparameters op de temperatuur. De huidige bipolaire radiofrequente ablatie wordt onderworpen aan een grondige analyse met als belangrijkste conclusie dat de warmtebron zeer lokaal is op spatiaal gebied. RFA kan in deze configuratie dus enkel een hulpmiddel zijn bij de behandeling van kleine tumoren. Verbeteringen zijn dus nodig en worden aangereikt in hoofdstuk 6. De bipolaire configuratie heeft immers het grote voordeel dat er twee ablatiezones tegelijkertijd ontstaan. Deze ablatiezones kunnen overlappen zodat er een temperatuursopbouw ontstaat in het tussenliggende gebied. Deze overlap blijkt in belangrijke mate be¨ınvloed te worden door het ge¨ınjecteerde stroomprofiel. We gebruiken een gepulseerde stroomvorm waarbij de invloed van de afschakeltijd geanalyseerd wordt. Deze afschakeltijd kunnen we veranderen door ofwel de aanschakeltemperatuur te variëren ofwel door te opteren voor een vaste afschakeltijd. Voor beide beslissingscriteria nemen we grote verschillen in efficiëntie waar tussen de verschillende afschakeltijden. Onze ontwikkelde numerieke methoden laten toe om voor elke patiënt of tumorgrootte na te gaan wat de optimale stroomvorm is voor een gegeven locatie van de naalden. Een patiëntspecifieke optimalisatie van bipolaire radiofrequente ablatie is dus in principe mogelijk.

7.2

Verder onderzoek

Een wijs gezegde luidt: “Hoe groter het eiland van kennis, hoe langer de kustlijn van verwondering.” We zijn ons bewust van het feit dat deze masterproef geen eindpunt vormt in de zin dat er nog op verschillende fronten mogelijkheden zijn tot verder onderzoek. We sluiten af met een summiere opsomming van enkele punten tot verder onderzoek. Vooreerst kan er vooruitgang geboekt worden op numeriek vlak: • Het wiskundig probleem dat opgelost diende te worden, werd gereduceerd tot een beginwaardenprobleem van eerste orde. Voor de oplossing hiervan, opteerden we in deze masterproef voor de modified midpoint-methode. Numeriek onderzoek kan uitgevoerd worden om te bepalen welke methode (Crank-Nicolson, Runge-Kutta, Bulirsch-Stoer, etc.) de gegeven differentiaalvergelijkingen het best oplost. • Verder fundamenteel wiskundig onderzoek kan uitgevoerd worden naar de stabiliteit van het gestelde probleem. Indien in de bio-warmtevergelijking ook convectietermen ge¨ıncludeerd worden om de invloed van grote bloedvaten te benaderen, is het noodzakelijk om meer extensieve analyses uit te voeren. Dit is zeker nodig indien het oplossingsdomein bestaat uit omgevingen waarin de convectie dominant is, terwijl de diffusietermen overheersen in een andere omgeving.

Hoofdstuk 7. Conclusie

95

• Het wiskundig model in dit proefschrift betreft een gekoppeld probleem. Een sensitiviteitsanalyse kan het belang aangeven van de koppeling tussen de systemen gekarakteriseerd door beide partiële differentiaalvergelijkingen. • Voor het oplossen van de Poisson-vergelijking werd een eindige differentiemethode toegepast. Een vergelijkende studie van de nauwkeurigheid van deze methode ten opzichte van deze van de eindige elementenmethode kan bijvoorbeeld uitgevoerd worden voor thermotherapie. • Ten slotte zouden we een stochastische oplossingsmethode kunnen opbouwen voor gegeven stochastische materiaalparameters. Fundamenteel onderzoek naar een dergelijke methode die de probabiliteitsdistributie van de temperatuursdistributie weergeeft kan van grote waarde zijn. Belangrijk hierbij is om een computationeel efficiënte methode op te bouwen die dit kan oplossen. Ten tweede kan er op biomedisch gebied, in termen van efficiëntie van thermotherapie, verder onderzoek uitgevoerd worden op volgende punten: • Andere naaldvormen zijn mogelijk waarbij gedacht kan worden aan naaldconfiguraties met spatiaal gedistribueerde stroominjecties die mogelijks tot een hogere temperatuur leiden in een groter spatiaal gebied. • In dit proefschrift werd een gepulste stroomvorm gebruikt voor de thermotherapie. Het is mogelijk om die stroomvorm iteratief te bepalen via een optimalisatie-algoritme. • Een meer homogene thermotherapie is mogelijk door het weefsel op te warmen via magnetische nanodeeltjes. Door het aanleggen van een tijdsafhankelijk magnetisch veld is het mogelijk om de deeltjes op te warmen ‘vanop afstand’. • Indien de tumor gelegen is nabij een groot bloedvat, zou de therapie slechter kunnen werken als gevolg van de convectie. Verder onderzoek kan erin bestaan om na te gaan hoe de therapieparameters aangepast moeten worden voor het behandelen van dergelijke tumoren. Ten laatste valt er ook op experimenteel gebied nog veel te onderzoeken: • Een uitvoerige experimentele validatie van de voorgestelde numerieke methoden kan nuttig zijn om de nauwkeurigheid van deze methoden te vergelijken met de werkelijkheid. Een mogelijkheid bestaat erin om dit onderzoek uit te voeren op varkenslevers (materiaalparameters zijn gelijkaardig aan die van de mens) en daarbij temperatuursmetingen uit te voeren door middel van thermokoppels. • Door de stroomvorm (‘on-line’) te controleren in werkelijke omstandigheden op basis van enerzijds de experimentele temperatuursmetingen en anderzijds op basis van het numeriek model, zouden we een verbetering van thermotherapie kunnen verkrijgen.

Bijlage A

De warmtevergelijking A.1

Afleiding van de warmtevergelijking

Beschouw een volume V van een bepaald materiaal dat een hoeveelheid warmte H (t) bevat [43]. Dan bestaat de volgende relatie tussen de temperatuur T (x, y, z, t) en de warmte H (t) in het volume V : ZZZ H (t) = ρ c T (x, y, z, t) d~v (A.1) V

h i kg en c de specifieke warmtecapaciteit van met ρ de massadichtheid van het materiaal m3 h i J het materiaal . De verandering in warmte is dan: kg K dH (t) = dt

ZZZ

ρ c Tt (x, y, z, t) d~v

(A.2)

V

Volgens de Wet van Fourier stroomt warmte van warme naar koude gebieden proportioneel met de temperatuursgradiënt [43]. Deze proportionaliteitsconstante is de thermische geleid baarheid k KWm . De warmte kan wegens het behoud van energie het volume V echter niet verlaten, tenzij door de rand van het volume ∂V . De verandering in warmte dH(t) is aldus dt ook gelijk aan: ZZ dH (t) = k (~n · ∇T (x, y, z, t)) dS (A.3) dt ∂V

Hierbij is ~n de naar buiten gerichte normaalvector aan het randoppervlak ∂V . Bijgevolg is: ZZZ ZZ ρ c Tt (x, y, z, t) d~v = k (~n · ∇T (x, y, z, t)) dS (A.4) V

∂V

Gebruik makend van de stelling van Gauss-Ostrogradsky, resulteert vgl. (A.4) in: ZZZ ZZZ ρ c Tt (x, y, z, t) d~v = ∇ · (k ∇T (x, y, z, t)) d~v V

(A.5)

V

In het volume V wordt vgl. (A.5) voorgesteld door de volgende partiële differentiaalvergelijking: ∂T (x, y, z, t) = ∇ · (k ∇T (x, y, z, t)) (A.6) ρc ∂t Deze partiële differentiaalvergelijking noemt men de warmtevergelijking [43]. 96

97

Bijlage A. De warmtevergelijking

A.2

Fundamentele oplossing van de warmtevergelijking in een homogeen medium

In een homogeen medium reduceert de warmtevergelijking (A.6) zich tot: ∂T k 2 = ∇ T = K ∇2 T ∂t ρc

(A.7)

h 2i Hierbij staat K gekend als de warmtediffusiviteit met als eenheid ms . De fundamentele oplossing van de eendimensionale warmtevergelijking in een homogeen medium wordt dus bekomen door het volgende beginwaardenprobleem op te lossen:   ∂u(x,t) − K ∇2 u (x, t) = 0 ∂t (A.8) u (x, t = 0) = δ (x)

De eenvoudigste manier om dit beginwaardenprobleem op te lossen is door een Fouriertransformatie toe te passen. Algemeen geldt dat de Fourier-transformatie van een absoluut integreerbare functie f (~x) met ~x = (x1 , · · · , xn ) gedefinieerd wordt als: Z Ff (~x) = fˆ ~k = f (~x) exp −i~k · ~x d~x (A.9) Rn

De inverse Fourier-transformatie definiëren we als: Z 1 fˆ ~k exp i~k · ~x d~k f (~x) = n (2π)

(A.10)

Rn

Een handige regel omtrent de Fourier-transformatie van de afgeleide van een absoluut integreerbare functie betreft: # " ∂ n f (~x) = (i kj )n Ff (~x) (A.11) F ∂xnj In beginwaardenprobleem (A.8) zijn we op zoek naar de fundamentele oplossing φ met als argument ~x = (x, t). De Fourier-getransformeerde van deze functie wordt genoteerd als: φˆ (k, ω) = Fφ (x, t)

(A.12)

We zien dus dat we vgl. (A.9) en (A.10) kunnen toepassen mits gebruik te maken van ~x = (x, t), ~k = (k, ω) en n = 2. Fourier-transformatie van beginwaardenprobleem (A.8) met φ (x, t) als fundamentele oplossing levert: ∂φ (x, t) (A.11) 2 − K ∇ φ (x, t) = δ (x) δ (t) = (i ω) Fφ (x, t) − K (i k)2 Fφ (x, t) = 1 F ∂t (A.12) = φˆ (k, ω) K k 2 + i ω = 1 (A.13) De Fourier-getransformeerde φˆ (k, ω) van de fundamentele oplossing bedraagt dus: φˆ (k, ω) =

1 K k2 + i ω

(A.14)

98


De fundamentele oplossing φ (x, t) kunnen we dan bekomen door op betrekking (A.14) een inverse Fourier-transformatie toe te passen: Z 1 1 (A.10) φ (x, t) = ei k x ei ω t dk dω 2 2 K k + iω (2π) =

We berekenen eerst

+∞ R

−∞

eiωt Kk2 +iω

1 (2π)2

R2 +∞ Z

dk e

+∞ Z

ikx

−∞

−∞

dω. Stellen we f (ω) =

ei ω t dω K k2 + i ω

eiωt , Kk2 +iω

(A.15)

dan heeft deze functie een

pool bij ω = iKk 2 in het bovenhalfvlak aangezien K > 0. Indien t > 0 dan kunnen we het lemma van Jordan toepassen op een halve cirkel in het bovenhalfvlak zodat: +∞ Z f (ω) dω = 2πi ResiKk2 f (ω)

−∞

e−K k = 2πi i

2

t

= 2πe−K k

2

t

(A.16)

Als t < 0, dan moeten we het lemma van Jordan toepassen op een halve cirkel in het +∞ R onderhalfvlak. Aangezien f (ω) geen polen heeft in het onderhalfvlak zal f (ω) dω = 0. −∞

Beide gevallen kunnen we samenvattend beschrijven door gebruik te maken van de Heavisidefunctie H (t): +∞ Z ei ω t 2 dω = 2πe−K k t H (t) (A.17) 2 K k + iω −∞

Substitutie van uitdrukking (A.17) in vgl. (A.15) levert: φ (x, t) =

=

=

H (t) 2π H (t) 2π

+∞ Z 2 ei k x e−K k t dk

−∞ +∞ Z −∞

"

2

exp i k x − K k t −

H (t) − x2 e 4K t 2π

ix √ 2 Kt

2

+

ix √ 2 Kt

" +∞ 2 # Z ix 2 √ exp i k x − K k t − dk 2 Kt

2 #

dk

(A.18)

−∞

√ √ Stellen we nu y = 2√i xKt − Ktk, dan is − Kt dk = dy. Verder is k = ±∞ ⇔ y = ∓∞. Met deze notatieverandering kunnen we vgl. (A.18) herschrijven als: φ (x, t) =

=

H (t) − x2 e 4K t 2π

−∞ Z exp −y 2

+∞

H (t) − x2 √ e 4K t 4π 2 Kt

1 √ dy − Kt

+∞ Z exp −y 2 dy

−∞

(A.19)

99


De gekende integraal van Poisson

√ exp −y 2 dy bedraagt π zodat:

+∞ R

−∞

φ (x, t) = √

√

π

x2

4π 2 Kt

e− 4 K t H (t)

(A.20)

We besluiten dat de fundamentele oplossing van de eendimensionale warmtevergelijking de volgende vorm aanneemt: x2 1 − H (t) (A.21) φ (x, t) = 1 exp 4Kt (4πKt) 2 Beschouwen we nu de omschrijving van de fundamentele oplossing van de n-dimensionale warmtevergelijking (~x = (x1 , . . . , xn )):  n P ∂u(~ x,t)  ∂2  u (~x, t) = 0  ∂t − K ∂xj 2 j=1 (A.22) n Q   δ (xi ) u(~x, t = 0) = δ (~x) = i=1

We bewijzen nu dat de volgende functie φ (~x, t) de oplossing is van (A.22): φ (~x, t) =

n Y

φ (xi , t)

(A.23)

i=1

Hierbij is φ (xi , t) =

1

1 (4πKt) 2

waardenprobleem:

x2i exp − 4Kt H (t) de oplossing van het 1-dimensionale begin  ∂φ(xi ,t) − K ∂ 2 φ(xi ,t) = 0 ∂t ∂xi 2 φ(x , t = 0) = δ (x ) i

(A.24)

i

Substitutie van vgl. (A.23) in de bovenste vergelijking van het beginwaardenprobleem (A.22) levert: n

X ∂2 ∂ φ (~x, t) − K φ (~x, t) = 0 ∂t ∂xj 2 j=1 ! ! n n n X Y ∂ Y ∂2 ⇐⇒ φ (xi , t) − K φ (xi , t) = 0 ∂t ∂xj 2 i=1 j=1 i=1     n n n n X Y X  ∂φ (xj , t) Y ∂ 2 φ (xj , t)   ⇐⇒ · φ (xi , t) − K   φ (xi , t) · =0 ∂t ∂xj 2 j=1

⇐⇒

n X j=1

i=1 i6=j

j=1

(A.25)

(A.26)

(A.27)

i=1 i6=j

  n ∂φ (xj , t) ∂ 2 φ (xj , t)  Y −K  φ (xi , t) · =0 ∂t ∂xj 2

(A.28)

i=1 i6=j

∂φ(x ,t)

∂ 2 φ(x ,t)

Aangezien ∂tj −K ∂xj 2j = 0 ∀j ∈ 1, . . . , n (vgl. (A.24)), zal φ (~x, t) inderdaad voldoen aan de eerste voorwaarde van beginwaardenprobleem (A.22). φ (~x, t) voldoet ook aan de

100

Bijlage A. De warmtevergelijking tweede voorwaarde van (A.22), immers: φ (~x, t = 0) = =

n Y

i=1 n Y

φ (xi , t = 0) δ (xi )

i=1

= δ (~x)

(A.29)

φ (~x, t) voldoet aan beide voorwaarden van beginwaardenprobleem (A.22) en is dus de oplossing van dit beginwaardenprobleem. In het bijzondere geval van drie dimensies (n = 3), reduceert de oplossing zich tot: 2 x + y2 + z2 1 H (t) (A.30) − φ (x, y, z, t) = 3 exp 4Kt (4πKt) 2

Bibliografie [1] 10voorBiologie - Hoofdstuk 38: //www.10voorbiologie.nl

Lever en nieren. [Online]. Available:

http:

[2] A. Bhattacharya and R. L. Mahajan, “Temperature dependence of thermal conductivity of biological tissues,” Physiological Measurement, vol. 24, no. 3, pp. 769 – 783, 2003. [3] J. Bischof, K. Christov, and B. Rubinsky, “A morphological study of cooling rate response in normal and neoplastic human liver tissue: Cryosurgical implications,” Cryobiology, vol. 30, no. 5, pp. 482 – 492, 1993. [4] M. Breen, M. Breen, K. Butts, L. Chen, G. Saidel, and D. Wilson, “MRI-guided Thermal Ablation Therapy: Model and Parameter Estimates to Predict Cell Death from MR Thermometry Images,” Annals of Biomedical Engineering, vol. 35, no. 8, pp. 1391 – 1403, 2007. [5] I. Chang, “Considerations for thermal injury analysis for RF ablation devices,” The Open Biomedical Engineering Journal, vol. 4, pp. 3 – 12, 2010. [6] I. Chang and U. Nguyen, “Thermal modeling of lesion growth with radiofrequency ablation devices,” BioMedical Engineering OnLine, vol. 3, no. 1, p. 27, 2004. [7] J. Chato, “Heat transfer to blood vessels,” Journal of Biomechanical Engineering, vol. 102, no. 2, pp. 110 – 118, 1980. [8] C.-C. R. Chen, M. I. Miga, and R. L. Galloway, “Optimizing electrode placement using finite-element models in radiofrequency ablation treatment planning,” IEEE Transactions on Biomedical Engineering, vol. 56, no. 2, pp. 237 – 245, 2009. [9] M. M. Chen and K. R. Holmes, “Microvascular contributions in tissue heat transfer,” Annals of the New York Academy of Sciences, vol. 335, pp. 137 – 150, 1980. [10] K. S. Cole and R. H. Cole, “Dispersion and absorption in dielectrics - I. Alternating current characteristics,” Journal of Chemical Physics, vol. 9, pp. 341 – 351, 1941. [11] S. Curley, “Radiofrequency ablation of malignant liver tumors,” The Oncologist, vol. 6, pp. 14 – 23, 2001. [12] G. Dahlquist and ˚ A. Bj¨ orck, Numerical Methods.

101

Prentice–Hall, Inc., 1974.

102

Bibliografie

[13] J. Davila, R. Morgan, Y. Shaib, K. McGlynn, and H. B. El-Serag, “Hepatitis C infection and the increasing incidence of hepatocellular carcinoma: a population-based study,” Gastroenterology, vol. 127, pp. 1372 – 1380, 2004. [14] T. de Baère, O. Risse, V. Kuoch, C. Dromain, C. Sengel, T. Smayra, M. G. E. Din, C. Letoublon, and D. Elias, “Adverse events during radiofrequency treatment of 582 hepatic tumors,” American Journal of Roentgenology, vol. 181, no. 3, pp. 695 – 700, 2003. [15] N. De Geeter, G. Crevecoeur, and L. Dupré, “An efficient 3-D eddy-current solver using an independent impedance method for transcranial magnetic stimulation,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 2, pp. 310 – 320, 2011. [16] H. De Schepper, Wiskundige Analyse II.

UGent, 2007 - 2008.

[17] G. D. Dodd, D. Napier, J. D. Schoolfield, and L. Hubbard, “Percutaneous radiofrequency ablation of hepatic tumors: Postablation syndrome,” American Journal of Roentgenology, vol. 185, no. 1, pp. 51 – 57, 2005. [18] I. dos Santos, D. Haemmerich, D. Schutt, A. F. da Rocha, and L. R. Menezes, “Probabilistic finite element analysis of radiofrequency liver ablation using the unscented transform,” Physics in Medicine and Biology, vol. 54, no. 3, pp. 627 – 639, 2009. [19] L. Dupré and G. Crevecoeur, Modelleren en simuleren van dynamische systemen. UGent, 2011 - 2012. [20] H. B. El-Serag, “Hepatocellular carcinoma,” New England Journal of Medicine, vol. 365, no. 12, pp. 1118 – 1127, 2011. [21] K. Friston, J. Ashburner, S. Kiebel, T. Nichols, and W. Penny, Eds., Statistical Parametric Mapping: The Analysis of Functional Brain Images. Academic Press, 2007. [22] C. Gabriel, S. Gabriel, and E. Corthout, “The dielectric properties of biological tissues: I. Literature survey,” Physics in Medicine and Biology, vol. 41, no. 11, pp. 2231 – 2249, 1996. [23] S. Gabriel, R. W. Lau, and C. Gabriel, “The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz,” Physics in Medicine and Biology, vol. 41, no. 11, pp. 2251 – 2269, 1996. [24] ——, “The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues,” Physics in Medicine and Biology, vol. 41, no. 11, pp. 2271 – 2293, 1996. [25] D. Haemmerich and D. Schutt, “RF ablation at low frequencies for targeted tumor heating: In vitro and computational modeling results,” IEEE Transactions on Biomedical Engineering, vol. 58, no. 2, pp. 404 – 410, 2011.

Bibliografie

103

[26] D. Haemmerich, S. T. Staelin, J. Z. Tsai, S. Tungjitkusolmun, D. M. Mahvi, and J. G. Webster, “In vivo electrical conductivity of hepatic tumours,” Physiological Measurement, vol. 24, no. 2, pp. 251 – 260, 2003. [27] D. Haemmerich and B. J. Wood, “Hepatic radiofrequency ablation at low frequencies preferentially heats tumour tissue,” International Journal of Hyperthermia, vol. 22, no. 7, pp. 563 – 574, 2006. [28] A. Jemal, F. Bray, M. M. Center, J. Ferlay, E. Ward, and D. Forman, “Global cancer statistics,” CA: A Cancer Journal for Clinicians, vol. 61, no. 2, pp. 69 – 90, 2011. [29] H. Kasugai, Y. Osaki, H. Oka, M. Kudo, and T. Seki, “Severe complications of radiofrequency ablation therapy for hepatocellular carcinoma: An analysis of 3,891 ablations in 2,614 patients,” Oncology, vol. 72 (Suppl. 1), pp. 72 – 75, 2007. [30] E. Kreyszig, Advanced Engineering Mathematics, 9th ed. 2006, ch. 12, pp. 588 – 590.

John Wiley & Sons, Inc.,

[31] D. S. K. Lu, S. S. Raman, P. Limanond, D. Aziz, J. Economou, R. Busuttil, and J. Sayre, “Influence of large peritumoral vessels on outcome of radiofrequency ablation of liver tumors,” Journal of vascular and interventional radiology, vol. 14, pp. 1267 – 1274, 2003. [32] D. S. K. Lu, N. C. Yu, S. S. Raman, P. Limanond, C. Lassman, K. Murray, M. J. Tong, R. G. Amado, and R. W. Busuttil, “Radiofrequency ablation of hepatocellular carcinoma: Treatment success as defined by histologic examination of the explanted liver,” Radiology, vol. 234, no. 3, pp. 954 – 960, 2005. [33] D. Miklavˇ ciˇ c, N. Pavˇ selj, and F. X. Hart, Wiley Encyclopedia of Biomedical Engineering. John Wiley & Sons, Inc., 2006, ch. Electric Properties of Tissues. [34] D. A. Nelson, “Invited editorial on “Pennes0 1948 paper revisited”,” Journal of Applied Physiology, vol. 85, no. 1, pp. 2 – 3, 1998. [35] Y. Ni, S. Mulier, Y. Miao, L. Michel, and G. Marchal, “A review of the general aspects of radiofrequency ablation,” Abdominal Imaging, vol. 30, pp. 381 – 400, 2005. [36] A. P. O’Rourke, D. Haemmerich, P. Prakash, M. C. Converse, D. M. Mahvi, and J. G. Webster, “Current status of liver tumor ablation devices,” Expert Review of Medical Devices, vol. 4, no. 4, pp. 523 – 537, 2007. [37] H. H. Pennes, “Analysis of tissue and arterial blood temperatures in the resting human forearm,” Journal of Applied Physiology, vol. 1, no. 2, pp. 93 – 122, 1948. [38] R. Pethig and D. B. Kell, “The passive electrical properties of biological systems: their significance in physiology, biophysics and biotechnology,” Physics in Medicine and Biology, vol. 32, no. 8, pp. 933 – 970, 1987. [39] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C++: The Art of Scientific Computing. Cambridge University Press, 2002.

104

Bibliografie

[40] S. Rossi, V. Ravetta, L. Rosa, G. Ghittoni, F. T. Viera, F. Garbagnati, E. M. Silini, P. Dionigi, F. Calliada, P. Quaretti, and C. Tinelli, “Repeated radiofrequency ablation for management of patients with cirrhosis with small hepatocellular carcinomas: A longterm cohort study,” Hepatology, vol. 53, no. 1, pp. 136 – 147, 2011. [41] H. Saleheen and K. Ng, “A new three-dimensional finite-difference bidomain formulation for inhomogeneous anisotropic cardiac tissues,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 1, pp. 15 – 25, 1998. [42] M. Slodiˇ cka and R. Van Keer, Numerieke Wiskunde.

UGent, 2009 - 2010.

[43] W. A. Strauss, Partial Differential Equations: An Introduction. Inc., 1992, ch. 1, p. 16.

John Wiley & Sons,

[44] R. Tateishi, S. Shiina, T. Teratani, S. Obi, S. Sato, Y. Koike, T. Fujishima, H. Yoshida, T. Kawabe, and M. Omata, “Percutaneous radiofrequency ablation for hepatocellular carcinoma. An analysis of 1000 cases,” Cancer, vol. 103, no. 6, pp. 1201 – 1209, 2005. [45] The Cancer of the Liver Italian Program Investigators, “A new prognostic system for hepatocellular carcinoma: A retrospective study of 435 patients,” Hepatology, vol. 28, no. 3, pp. 751 – 755, 1998. [46] J. W. Valvano, J. R. Cochran, and K. R. Diller, “Thermal conductivity and diffusivity of biomaterials measured with self-heated thermistors,” International Journal of Thermophysics, vol. 6, pp. 301 – 311, 1985. [47] N. Van den Bergh, Wiskundige Methoden.

UGent, 2010 - 2011.

[48] S. Weinbaum and L. M. Jiji, “A new simplified bioheat equation for the effect of blood flow on local average tissue temperature,” Journal of Biomechanical Engineering, vol. 107, no. 2, pp. 131 – 139, 1985. [49] S. Weinbaum, L. X. Xu, L. Zhu, and A. Ekpene, “A new fundamental bioheat equation for muscle tissue: Part I - Blood perfusion term.” Journal of Biomechanical Engineering, vol. 119, no. 3, pp. 278 – 288, 1997. [50] S. L. Wong, P. B. Mangu, M. A. Choti, T. S. Crocenzi, G. D. Dodd, G. S. Dorfman, C. Eng, Y. Fong, A. F. Giusti, D. Lu, T. A. Marsland, R. Michelson, G. J. Poston, D. Schrag, J. Seidenfeld, and A. B. Benson, “American Society of Clinical Oncology 2009 clinical evidence review on radiofrequency ablation of hepatic metastases from colorectal cancer,” Journal of Clinical Oncology, vol. 28, no. 3, pp. 493 – 508, 2010. [51] X.-Y. Yin, X.-Y. Xie, M.-D. Lu, H.-X. Xu, Z.-F. Xu, M. Kuang, G.-J. Liu, J.-Y. Liang, and W. Y. Lau, “Percutaneous thermal ablation of medium and large hepatocellular carcinoma,” Cancer, vol. 115, no. 9, pp. 1914 – 1923, 2009. [52] L. Zhu, Standard Handbook of Biomedical Engineering and Design. McGraw-Hill, 2003, ch. 2, pp. 2.1 – 2.29.

Lijst van tabellen 2.1 2.2 2.3 2.4

Materiaaleigenschappen van levercellen bij toepassing van RFA [27]. . . . . . Gemiddeld aantal elektrode-inserties in functie van de tumorgrootte [44]. . . . Verband tussen tumornecrose en de aanwezigheid van een groot bloedvat [32]. Verband tussen tumornecrose en de tumorgrootte [51]. . . . . . . . . . . . . .

3.1

Materiaaleigenschappen van de weefsels betrokken in de numerieke simulatie van RFA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fysische constanten die de bloedperfusieterm karakteriseren [8]. . . . . . . . . Overzicht van de parameters gebruikt in het Arrheniusmodel voor thermische schade van biologische weefsels. . . . . . . . . . . . . . . . . . . . . . . . . . . Definitie van de statistische parameters die de uitkomst van een RFA-procedure karakteriseren. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2 3.3 3.4

5.1 5.2 5.3

Materiaaleigenschappen van de verschillende weefsels die aanwezig zijn in de simulatieomgeving. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Afmetingen die de verschillende geometrieën karakteriseren om de invloed van discretisatie na te gaan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Benodigde rekentijd voor de simulatie van de verschillende geometrieën. . . .

105

19 21 22 22

26 26 42 44

69 75 75

Lijst van figuren 2.1 2.2 2.3 2.4 2.5 2.6 2.7 3.1 3.2

Schematische voorstelling van een elementair volume weefsel dat slechts één slagader-ader paar bevat. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Schematische opbouw van een leverlobje [1]. . . . . . . . . . . . . . . . . . . . Diëlektrisch spectrum van biologische weefsels. . . . . . . . . . . . . . . . . . Elektrische geleidbaarheid van normale levercellen en tumorcellen in functie van de frequentie [26]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Experimentele verwarmings- en afkoelingscyclus tot 90 ◦ C waarbij de relatie tussen k en T niet-lineaire effecten vertoont [2]. . . . . . . . . . . . . . . . . . Schematische voorstelling van de opstelling gebruikt in huidige RFA-procedures. Schematische voorstelling van een nieuwe configuratie voor RFA. . . . . . . . Dwarsdoorsnede van een simulatieomgeving. . . . . . . . . . . . . . . . . . . . Overzicht van de globale werking van het numeriek model voor radiofrequente ablatie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Schematische voorstelling van de analytische situatie aangaande de bolkapelektrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Grafische weergave van het stroomdichtheidsprofiel ge¨ınjecteerd door een bolkapelektrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Numerieke resultaten van een elektromagnetisch deelprobleem. . . . . . . . . 4.4 Vergelijking van de genormeerde numerieke en analytische oplossingen in het geval van de delta-bron. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Logaritme van het relatieve verschil tussen de genormeerde numerieke en analytische oplossingen in het geval van de delta-bron. . . . . . . . . . . . . . . . 4.6 Vergelijking van de genormeerde numerieke en analytische oplossingen in het geval van de bolkapelektrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Logaritme van het relatieve verschil tussen de genormeerde numerieke en analytische oplossingen in het geval van de bolkapelektrode. . . . . . . . . . . . . 4.8 Grafische weergave van het functievoorschrift voor u (xi ). . . . . . . . . . . . 4.9 Numerieke en analytische oplossing in het punt (26, 26, 26) binnen de centrale regio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Verhouding van de numerieke en analytische oplossing in het punt (26, 26, 26) binnen de centrale regio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 12 16 17 19 20 23 26 47

4.1

106

52 53 56 57 57 58 58 60 63 63

Lijst van figuren 4.11 Numerieke en analytische oplossing in het punt (11, 11, 11) buiten de centrale regio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Verhouding van de numerieke en analytische oplossing in het punt (11, 11, 11) buiten de centrale regio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Numerieke en analytische oplossing op het tijdstip t = 100 s in de snede {y = 0, z = 0}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Verhouding van de numerieke en analytische oplossing op het tijdstip t = 100 s in de snede {y = 0, z = 0}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Numerieke en analytische oplossing in het punt (26, 26, 26) voor verschillende stapgrootten in de tijd. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Numerieke oplossing in de snede z = −0.025 m waarbij de evolutie van de numerieke instabiliteit wordt gevisualiseerd. . . . . . . . . . . . . . . . . . . . 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13

6.1 6.2 6.3 6.4

Dwarsdoorsnede van de simulatieomgeving. . . . . . . . . . . . . . . . . . . . Driedimensionale impressie van de simulatieomgeving. . . . . . . . . . . . . . Vergelijking van het temperatuursverloop in het punt (51, 45, 51) in aan- of afwezigheid van een tumor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Logaritme van het verschil in potentiaal tussen beide configuraties in de snede x = 51. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Vergelijking van het temperatuursverloop in de snede {x = 51, z = 51} op het tijdstip 1.9 s in aan- of afwezigheid van een tumor. . . . . . . . . . . . . . . . Vergelijking van het temperatuursverloop in het punt (51, 45, 51) voor σ1 = 0 K−1 en σ1 = 0.0015 K−1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spatiale weergave van het verschil in temperatuur op de tijdstippen van maximale temperatuur voor σ1 = 0 K−1 en σ1 = 0.0015 K−1 . . . . . . . . . . . . . Norm van het temperatuursverschil binnen eenzelfde z-snede in functie van de co¨ ordinaat in de z-richting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperatuursverloop in de y-richting langsheen de naalden op het tijdstip van maximale temperatuur voor de verschillende discretisatieniveaus. . . . . . . . Vergelijking van het temperatuursverloop in het punt (51, 45, 51) met of zonder bloedperfusie. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Norm van de temperatuursverschillen in de snede zf = 53 in functie van de tijd. Vergelijking van het temperatuursverloop in het punt (51, 45, 51) met of zonder modulatie van de Pennes bio-warmtevergelijking. . . . . . . . . . . . . . . . . Vergelijking van de overlevingsgraad in het punt (51, 45, 51) met of zonder modulatie van de Pennes bio-warmtevergelijking. . . . . . . . . . . . . . . . . Referentiestroomvorm corresponderend met het beslissingscriterium inzake de temperatuur. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Stroomprofiel gedurende een RFA-procedure van 10 pulsen met aanduiding van de pulslabels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Overlevingsgraad op locaties dichtbij de naald en het centrum in functie van de pulslabels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Temperatuur na elke puls ter hoogte van het centrum (26, 26, 26). . . . . . . .

107

64 64 65 65 66 67 68 69 71 71 72 73 74 74 75 76 77 78 78

81 82 82 83

Lijst van figuren 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17

Temperatuursverval in de tijd indien geen stroom geleverd wordt. . . . . . . . Temperatuur na elke puls ter hoogte van het centrum voor verschillende aanschakeltemperaturen. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Overlevingsgraad ter hoogte van het centrum bij elk pulslabel en voor verschillende aanschakeltemperaturen. . . . . . . . . . . . . . . . . . . . . . . . . . . Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende aanschakeltemperaturen. . . . . . . . . . . . Overlevingsgraad ter hoogte van het centrum bij elk pulslabel en voor verschillende afschakeltijden. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spatiale weergave van de overlevingsgraad in de snede z = 26. . . . . . . . . . Temperatuur ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden. . . . . . . . . . . . . . . . . . . . Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden. . . . . . . . . . . . . . . . . Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede x = 26. Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede y = 26. . Spatiale weergave van de overlevingsgraad na 500 pulsen in de snede z = 26. . Overlevingsgraad ter hoogte van het centrum in functie van het aantal seconden stroom en voor verschillende afschakeltijden. . . . . . . . . . . . . . . . .

108 84 85 85 86 87 87 88 89 89 90 91 91 92

Ontwikkeling van numerieke methoden voor thermotherapie van leverkanker

Recommend Documents