RESEARCH ARTICLE
Fast Estimation of the Vascular Cooling in RFA Based on Numerical Simulation
T. Kröger^{*, 1}, T. Pätz^{1}, I. Altrogge^{2}, A. Schenk^{1}, K.S. Lehmann^{3}, B.B. Frericks^{3}, J.-P. Ritz^{3}, H.-O. Peitgen^{1, 2}, T. Preusser^{1, 4}
Article Information
Identifiers and Pagination:
Year: 2010Volume: 4
First Page: 16
Last Page: 26
Publisher Id: TOBEJ-4-16
DOI: 10.2174/1874120701004010016
Article History:
Received Date: 17/9/2009Revision Received Date: 22/11/2009
Acceptance Date: 28/12/2009
Electronic publication date: 4/2/2010
Collection year: 2010
open-access license: This is an open access article licensed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/) which permits unrestricted, non-commercial use, distribution and reproduction in any medium, provided the work is properly cited.
Abstract
We present a novel technique to predict the outcome of an RF ablation, including the vascular cooling effect. The main idea is to separate the problem into a patient independent part, which has to be performed only once for every applicator model and generator setting, and a patient dependent part, which can be performed very fast. The patient independent part fills a look-up table of the cooling effects of blood vessels, depending on the vessel radius and the distance of the RF applicator from the vessel, using a numerical simulation of the ablation process. The patient dependent part, on the other hand, only consists of a number of table look-up processes. The paper presents this main idea, along with the required steps for its implementation. First results of the computation and the related ex-vivo evaluation are presented and discussed. The paper concludes with future extensions and improvements of the approach.
1. INTRODUCTION
Cancer is currently a leading cause of death [1]. Over 50 % of all patients suffering from colon cancer develop hepatic metastases [2]. Since surgical resection is often not applicable due to size, number, and location of the tumors, or due to the general constitution of the patient, minimally invasive therapies, such as radiofrequency (RF) ablation, have become very popular in the last years for both primary hepatic tumors and hepatic metastases [3].
In RF ablation, a needle shaped probe is inserted into the tumor. The probe is then connected to an electric generator, which induces an electric current into the tissue. By Joule’s law, heat is produced, and above a certain critical temperature cell proteins denaturate, causing cell death.
Intervention planning is hampered by the fact that blood vessels in the vicinity of the RF probe can lead to a cooling effect and thus decrease the coagulation size [4]. Next to the area where the tissue is completely destroyed, there is a transition zone in which the tissue is still heated up, but the heat is not sufficient to cause cell death. There are some indications that possibly remaining tumor cells in this transition zone might become more aggressive and resistant towards further treatment than they would have been if no RF ablation had been performed: It is known that increased temperatures give rise to expression of heat shock proteins [5, 6], which may in turn later-on protect the cells against heat-induced apoptosis [7]. Also, the transition zone develops an increased activity of matrix metalloproteinase [8], which is known to be connected to an increased risk of tumor recurrence in colorectal cancer [9]. Although the clinical relevance of these correlations is yet unknown and there are also a number of indications for the converse point of view (see also the overview paper of Mulier et al. [10]), we argue that incomplete tumor ablation should be avoided whenever possible.
From the slices of a CT scan, it is typically even hard to tell how distant a tumor is located from a large blood vessel, not to mention the cooling effect of this blood vessel on the necrosis. This causes a need for computer assisted planning tools and simulations of the ablation, including the vascular cooling effect. Such simulations are developed by various authors. The way in which vascular cooling is taken into account varies considerably: Simple models use heuristics to quantify the cooling effects [11] or do not consider patient individual vascular systems at all [12]. More advanced approaches model individual vessels by a fixed heat sink [13, 14], and some authors even take blood flow together with heat exchange between vessel and tissue into account [15-17]. For a comprehensive review, see the overview article by Berjano [18]. However, the more precise cooling is considered, the slower the computations are, and, in fact, all models—except those that are based on heuristics—are far too slow to be usable in clinical practice.
In the current publication, we present a novel approach to overcome this dilemma. The idea is to separate the computationally intensive part of the simulation from the patient dependent part as much as possible. Computations that are not dependent on the individual anatomy of the patient can be performed in advance and, if necessary, on high-performance computers. For a given RF applicator and generator type, our idea is to precompute the cooling effect of blood vessels, depending on the vessel radius and its distance to the applicator, and to store these results in a table. Once the table has been generated, the cooling effects for any individual case can be estimated and visualized in real time. Our method points out tumor regions that would be underablated due to vascular cooling and can thus help the attending physician to choose a suitable applicator placement or to decide whether additional steps like a Pringle manoeuvre or chemoembolization should be applied.
2. MATERIALS AND METHODS
In Sect. 2.1, we describe the basic concept of the method, consisting of a quantitative description of necrosis zones in the presence of a vessel as well as an extension model for the case of a complete vascular tree. Section 2.2 briefly addresses the creation of an adaptive look-up table. The mathematical model of RF ablation that we use in this work is presented in Sect. 2.3. along with its numerical approximation in Sect. 2.4.. The setup for our example computations and a first ex-vivo evaluation are described in Sect. 2.5.
2.1. Quantitative Description of Necrosis
2.1.1. Single Vessel Model
Let us first simplify the general setting by considering the situation of one straight blood vessel of infinite length, oriented parallel to the RF applicator. The cooling effect of this vessel—for a fixed applicator type and generator setting—is assumed to depend on the radius r_{ves} of the vessel and on the distance d_{ap} of the applicator from the vessel. (Other influence factors, such as the vessel type, will be respected in our forthcoming research.) We quantify the cooling effect by stating the size r_{nec} of the necrosis, measured from the applicator’s center, in all directions orthogonal to the applicator axis.
Using a forward simulation (cf. Sect. 2.3), we establish a connection between the vessel radius r_{ves}, the distance d_{ap} of the applicator from the vessel, and the resulting necrosis size r_{nec}, which is a function of the direction, expressed as the angle α_{nec} between the direction under consideration and the direction in which the vessel center is located. For reasons of symmetry, it suffices to consider α_{nec}∈[0,π]. We formally denote this relation as a function
$:\left(0,\left(\right)\left(0,\left(\right)\right)\left(0,\right)\right)\left(0,\left(\right)\left({r}_{\mathrm{ves}},{d}_{\mathrm{ap}},{}_{\mathrm{nec}}\right){r}_{\mathrm{nec}}\right)$
see also Fig. (1). For a precise definition we refer the reader to Sect. 2.3.
Fig. (1). Schematic sketch of the quantities that are related by the function Ф. |
Fig. (2). Schematic sketch of the situation in the case of an encasement. |
Fig. (3). Location of the points a_{i} and b_{i} (for some value of i) as well as a^{*} and b^{*}. |
Fig. (5). In the absence of near-by vessels, the estimated necrosis consists of the union of balls around all the points on the electrode. |
Fig. (6). Schematic sketch of the geometric setup. |
Fig. (7). Plot of the necrosis (solid line) for selected values of the vessel radius r_{ves} and the applicator distance d_{ap} , as reconstructed from the look-up table, i. e. from the function $\stackrel{}{}$ . For reasons of symmetry, we only show the top half plane. For clarity, the corresponding vessel (gray) is also shown. Given quantities are in millimeters. An encasement occured in the top left case—the (very small) tube of unaffected tissue (thickness θ_{tube} = 0.05mm ) is not shown. Framed plots correspond to the ex-vivo experiments shown in Figs. (10-12). |
Fig. (11). Result of an ex-vivo experiment with r_{ves} = 2.5mm and d_{ap} = 10mm . |
Fig. (12). Result of an ex-vivo experiment with r_{ves} = 5mm and d_{ap} = 5mm . |
Fig. (13). Schematic of different cooling effects induced by portal venous (left) and hepatic venous (right) vessel systems (qualitatively). |
In the case of an encasement of the vessel, we measure the distance to the outer-most boundary of the necrosis area. This measurement ignores the fact that there will be a thin layer of unaffected tissue around the vessel. To compensate for this, we additionally measure the thickness θ_{tube} of this tube and hence employ another function
$:\left(0,\left(\right)\left(0,\left(\right)\left(0,\right)\right)\right)\left({r}_{\mathrm{ves}},{d}_{\mathrm{ap}}\right){}_{\mathrm{tube}}$
see Fig. (2), where we approximate the tube by a circular cylinder.
2.1.2. Full Model
Now, let the complete vascular system be given as a set of straight vessel segments with end points a_{i}, b_{i} ∈ R^{3} and radii r_{i} > 0 where i=1,...,N (such a charaterization can be obtained using the methods presented earlier [19, 20]), and let a^{*} and b^{*} denote the applicator’s active zone’s end points and r^{*} the applicator’s radius (cf. Fig. 3). We make the simplification that the cooling effect of the blood vessels is additive in the sense that the necrosis achieved in the presence of two blood vessels A and B is given by the intersection of those achieved in the presence of only vessel A and only vessel B – see Fig. (4). For clinically relevant settings, this can be assumed to be a sufficiently good approximation.
In the following, we denote by [a,b] the straight line from a ∈ R^{3} to b ∈ R^{3} and by ∠abc the angle between the lines $\stackrel{}{\mathrm{ab}}$ and $\stackrel{}{\mathrm{bc}}$ . Additionally, by d_{r,s} (a,b) we denote the distance of a ball of radius r around a from a ball of radius s around b, i. e.
${d}_{r,s}\left(a,b\right)=|ab|rs$
Let ξ ∈ R^{3} be a given point in the tissue, and we want to decide whether the tissue at this point will be destroyed by the ablation or not. For each choice of a point x ∈ [a^{*} ,b^{*}] on the applicator’s active zone and a point y ∈ [a_{i} ,b_{i}] in the i’th vessel segment, we consider the plane containing ξ, x, and y and apply the considerations described in Sect. 2.1.1 in this plane.
The transformation to that plane is not actually computed since the distances of the points ξ, x, and y from each other together with the angle in x suffice to obtain the required information. ξ is destroyed if it is contained in the necrosis described by the function Ф, but not in the tube described by the function Ψ; that is, if
$\left(x\right)\left(r,{d}_{{r,r}^{}}\left(y,x\right),\mathrm{yx}\right)\mathrm{and}\left(y\right)\left(r,{d}_{{r,r}^{}}\left(y,x\right)\right)$
We therefore define
(1a) ${C}_{y,{r,r}^{}}^{\left(1\right)}\left(x\right)=\left({R}^{3}:\left(x\right)\left(r,{d}_{{r,r}^{}}\left(y,x\right),\mathrm{yx}\right)\right)$
(1b) ${C}_{{y,r,r}^{}}^{\left(2\right)}\left(x\right)=\left({R}^{3}:\left(y\right)\left(r,{d}_{{r,r}^{}}\left(y,x\right)\right)\right)$
so that the set ${C}_{{y,r,r}^{}}^{\left(1\right)}\left(x\right){C}_{{y,r,r}^{}}^{\left(2\right)}\left(x\right)$ describes the necrosis for fixed x and y.
Now, let y vary over the vascular tree. This corresponds to taking several vascular segements into account and hence reduces the necrosis to the intersection of the corresponding sets. Finally, let x vary along the active zone’s centerline. Since longer electrodes generate longer coagulation zones, we take the union of all the respective necrosis sets. The final formula for the necrosis hence reads
(1c) $C=\underset{x\left({a}^{},{b}^{}\right)}{}\underset{i=1}{\overset{N}{}}\underset{y\left({a}_{i},{b}_{i}\right)}{}\left({C}_{y,{r}_{i},{r}^{}}^{\left(1\right)}\left(x\right){C}_{y,{r}_{i},{r}^{}}^{\left(2\right)}\left(x\right)\right)$
For the case where all vessel segments are too distant from the applicator to have any effect, this results in the union of balls with equal radius around all the points on the active zone, see Fig. (5).
2.2. Adaptive Creation of Look-Up Table
Let us assume that the functions Ф and Ψ (cf. Sect. 2.1) are known in the sense that for given r_{ves}, d_{ap} and α_{nec} , the values Ф(r_{ves}, d_{ap}, α_{nec}) and Ψ(r_{ves}, d_{ap}) can be computed, but this computation is too slow for clinical application. The idea is to approximate Ф and Ψ by some functions $\stackrel{}{}$ and $\stackrel{}{}$ , respectively, with the following properties:
- The construction of $\stackrel{}{}$ and $\stackrel{}{}$ (for a given tolerance) requires preferably few evaluations of Ф and Ψ.
- Once $\stackrel{}{}$ and $\stackrel{}{}$ have been computed, they can be evaluated very fast.
An easy possiblity consists of the evaluation of the functions on a uniform grid and storage of the function values in an array, together with a piecewise multilinear interpolation (that is, trilinear for Ф and bilinear for Ψ). However, it can be expected that this will result in many unnecessary evaluations since the functions might be approximately multilinear in some regions and more essentially non-multilinear in others.
For that reason, we use an adaptive hierarchic grid technique as described by Bungartz and Dirnstorfer [21], consisting of a hierachical ansatz space of piecewise multilinear functions. This produces a fine grid (i. e. small pieces on which $\stackrel{}{}$ and $\stackrel{}{}$ are multilinear) in regions where this is necessary only. However, in the α_{nec} direction, we use a full grid, since all values can be obtained from the same simulation result here.
Analogously to (1), we define approximations $\stackrel{}{C},{\stackrel{}{C}}_{y,r,{r}^{}}^{\left(1\right)},\text{and}{\stackrel{}{C}}_{y,r,{r}^{}}^{\left(2\right)}.$ . We emphasize that no special properties of any particular forward simulation model are utilized here.
2.3. A Mathematical Model of RF Ablation
In this section, we briefly describe the forward simulation model that we used to produce the results shown in Sect. 3, i. e. the employed functions Ф and Ψ (in the notation of Sect. 2.1).
As mentioned briefly in the introduction, a broad variety of simulation models for RF ablation exist, and they differ not only in the way in which they consider vascular cooling. They range from simple geometric models that use ellipsoids to describe the necrosis [11] over time dependent PDE models with constant coefficients [13, 16] to complex models that take cell water evaporation and temperature dependent material parameters into account [14, 15, 17]. Some authors exploit rotational symmetry around the applicator axis [12], which is of course impossible if the aim is to consider patient individual vascular systems.
Since the precomputation typically involves several hundred of calls of the forward simulation, we currently use a simplified forward simulation model here. It assumes constant material parameters and is based on a steady state solution. Since also the electric conductivity is assumed to be equal in the parenchyma and the blood vessel, the computation of the electric power density is independent of the location and the geometry of the vessel. Hence, for a given applicator type and generator setting, the potential equation (see (2) below) has to be solved only once, whereas the steady state bio-heat equation (see (4) below) has to be solved for varying blood vessel configurations. We further simplify our model by assuming the applicator to be placed in a direction parallel to the blood vessel and then computing the steady state bio-heat transfer equation in two space dimensions only. Note that for reasons that will become clear below, the potential and power computation must, in contrast, be computed in three space dimensions.
More specific, we consider a cuboid shaped compuational domain Ω ⊂ R^{3} that is chosen large enough such that any influences from outside can be neglected. Without loss of generality, we assume the centerlines of both the applicator and the blood vessel to be aligned in x_{3} direction and have an x_{2} coordinate of 0, cf. Fig. (6). The center of the applicator’s active zone is placed in the origin; the x_{1} position of the blood vessel is chosen to be negative and such that the distance between blood vessel and applicator matches the prescribed value d_{ap}. The open subdomains of Ω occupied by blood vessel, electrodes, and applicator are denoted by Ω_{ves}, Ω_{±}, and Ω_{pr} . Here, one of the electrodes is arbitrarily chosen to be Ω_{+} and the other Ω_{-} , both of which are subsets of Ω_{pr} ; in the case of a monopolar probe, only Ω_{+} is used. Further, we let
$\begin{array}{cc}{}_{\mathrm{out}}=& {}_{\mathrm{ves}}={}_{\mathrm{ves}}\\ {}_{}={}_{}{}_{\mathrm{pr}}& {}_{\mathrm{iso}}=\left({}_{\mathrm{pr}}\right){}_{}\end{array}$
Since the electric potential φ(x) can be considered as quasistatic [12], we can model it as the solution of the elliptic boundary value problem [14]
(2d) $n\left(x\right)\left(x\right)=\frac{n\left(x\right)\left(sx\right)}{{\left(sx\right)}^{2}}\left(x\right)\mathrm{on}{}_{\text{out}}$
Here, σ>0 is the electric conductivity, s is the center of the applicator’s active zone, and n(x) is the outward normal vector. It is a well-known fact that this problem has a unique weak solution ${H}^{1}\left({\stackrel{}{}}_{}\right)$ where H^{1} = H^{1,2} denotes the Sobolev space of functions whose weak first derivatives are square integrable [22].
The electric power density p is now given by [14]
$p\left(x\right)=\frac{{p}_{\mathrm{eff}}}{{p}_{\mathrm{total}}}{\left(\left(x\right)\right)}^{2}\mathrm{where}$
${p}_{\text{total}}=\underset{\stackrel{}{{}_{\text{pr}}}}{}{\left(\left(x\right)\right)}^{2}\mathrm{dx}\left(\mathrm{normalization\; factor}\right)$
${p}_{\text{eff}}=\frac{{4p}_{setup}{R}_{\text{tis}}{R}_{I}}{{\left({R}_{\text{tis}}{R}_{\text{I}}\right)}^{2}}\text{}\left(\text{effective power}\right)$
${R}_{\mathrm{tis}}=\frac{{U}^{2}}{{p}_{\mathrm{total}}}\text{}\left(\text{tissue impedance}\right)$
$\text{U}=\left(\begin{array}{c}1V\text{monopolar}\\ 2V\text{bipolar}\end{array}\right)\text{}\left(\text{difrence of the values of}\text{on electrodes}\right)$
and the remaining quantities are properties and settings of the electric generator, namely its inner resistance R_{I} and the setup power p_{setup}. (These formulae arise from an equivalent circuit that models the generator by a constant voltage source and an inner resistance. Even though we use a steady state model, it is important to capture this power scaling since it is a fact that the effective power is generally lower than the setup power, cf. Sects. 2.5.2 and 3.2).
It is well-known [23] that $p{L}^{1}\left(\stackrel{}{}\mathrm{pr}\right){H}^{1}\left(\stackrel{}{}\mathrm{pr}\right)$ , and if we extend p to Ω by setting p|_{Ωpr} =0, we have $p{L}^{1}\left(\right){H}^{1}\left(\right).$
Performing these computations in two (rather than in three) space dimensions would result in a mismatch of units: p_{total} would be an integral over a two-dimensional domain only, hence have a unit that differs by a factor of a meter from the correct unit. In consequence, this would R_{tis} cause to have a wrong unit as well, and finally the equation for p_{eff} would not be well-defined any more, since two quantities with different units would be added in the denominator. The physics of electric fields and electric potentials behaves differently in two dimensions than in three.
The temperature equation, on the other hand, can easily be performed two-dimensional. To do this, we let $\stackrel{}{}$ denote the intersection of Ω and the x_{1}-x_{2} plane as a subset of R^{2}, i. e.$\stackrel{}{}=\left(\left({x}_{1},{x}_{2}\right):\left({x}_{1},{x}_{2},0\right)\right)$ (see again Fig. 6) and analogously for Ω_{ves}, Ω_{pr}, and Г_{out}. Further, we let $\stackrel{}{p}{H}^{1}\left(\stackrel{}{}\right)$ be given by $\stackrel{}{p}=P\left(p\right)$ where
is a suitable projection to be specified below. The remaining modelling is completely performed in two space dimensions, but for simplicity of notation, we will omit the ^{~} symbol from now on.
The temperature distribution T(x) is modelled by the two-dimensional steady state bio-heat equation [14, 24, 25]
(4a) $T\left(x\right)=p\left(x\right)v\left({T}_{\mathrm{body}}T\left(x\right)\right)\mathrm{in}\text{}\left({\stackrel{}{}}_{\mathrm{ves}}\right)\left({\stackrel{}{}}_{\mathrm{pr}}\right)$
(4c) $T\left(x\right)={T}_{\mathrm{body}}\text{on}{\stackrel{}{}}_{\mathrm{ves}}{\stackrel{}{}}_{\mathrm{pr}}$
where λ>0 is the thermal conductivity, T_{body} is the native body temperature, and ν quantifies the cooling due to small vessels, which are assumed to pervade the complete tissue. This is again an elliptic boundary value problem with unique weak solution $T{H}^{1}\left(\right)$ .
Once T is known, we let Ω_{coag} ⊂ Ω denote the set of coagulated tissue,
${}_{\mathrm{coag}}=\left(x:T\left(x\right){T}_{\mathrm{crit}}\right)$
(where T_{crit} is a critical temperature above which the tissue is assumed to be destroyed).
Finally, we define the necrosis size r_{nec} in a direction α_{nec} and the thickness θ_{tube} of the uneffected tissue tube according to the size and shape of Ω_{coag} as described in Sect. 2.1.1. Since, however, a point evaluation of a function in H^{1} is in general not defined, the set Ω_{coag} is only defined up to sets of measure zero, so that the precise definition of r_{nec} and θ_{tube} is given by
${r}_{\mathrm{nec}}=\mathrm{sup}\left(r0:{}_{0}\left({B}_{}\left({\mathrm{re}}_{{}_{\mathrm{nec}}}\right){}_{\mathrm{coag}}\right)0\right)$
${}_{\mathrm{tube}}=\mathrm{sup}\left(r0:\left({B}_{r{r}_{\mathrm{ves}}}\left({x}_{\mathrm{ves}}\right){}_{\mathrm{coag}}\right)=0\right)$
where e_{α} = (cos α, sin α)and x_{ves} is the position of the vessel’s centerline in Ω (i. e.x_{ves} = (–r_{ves}–d_{ap},0)), B_{ε}(x) denotes a ball of radius ε around x, and $||$ denotes the Lebesgue measure of a set.
Now that the description of our functions Ф and Ψ is finished, we would like to add some remarks:
- The described forward simulation model provides only a rough approximation of the ablation outcome. This is, however, not a substantial drawback since any other model can be plugged in easily, thus any level of precision can be achieved without slowing down the patient individual estimation.
- It is well-known that a chemical process like the denaturation of liver cells is not only a function of temperature, but in fact also depends on the duration for which a given temperature is applied. Hence, a better model for the set Ω_{coag} would be to apply the Arrhenius formalism [12, 14, 26]. However, a steady-state solution does obviously not supply the required time information.
- A possible choice of the projection P (cf. (8)) is
$\left({P}_{}p\right)\left({x}_{1},{x}_{2}\right)=\frac{\underset{}{\overset{}{}}p\left({x}_{1},{x}_{2},{x}_{3}\right){}_{}\left({x}_{3}\right){\mathrm{dx}}_{3}}{\underset{}{\overset{}{}}{}_{}\left({x}_{3}\right){\mathrm{dx}}_{3}}$
where 2ℓ is the length of the applicator’s active zone and ${}_{}{C}^{}\left(R\right)$ is a non-negative function with compact support in [-ℓ,ℓ] and the property that ${}_{}{}_{\left(,\right)}1$ . It is straight forward to prove that P_{ε} is well-defined, but the norm does not remain bounded for ε→0. In other words,
$\left({P}_{0}p\right)\left({x}_{1},{x}_{2}\right)=\frac{1}{2}\underset{}{\overset{}{}}p\left({x}_{1},{x}_{2},{x}_{3}\right){\mathrm{dx}}_{3}$
does apparently not provide a mapping ${H}^{1}\left(\right){H}^{1}\left(\stackrel{}{}\right)$ , and although P_{ε} does provide such a mapping for any ε>0, the result might theoretically depend essentially on the (arbitrary) choice of ε. In the view of remark number 1, we argue that for a full (three-dimensional) forward model, no such mapping will be required any more anyway. Also, the numerical approximation p_{h} of the power density p (see Sect. 2.4 below) is an element of L^{∞}(Ω), so that P_{0}(p_{h}) is well-defined, and so is even P_{*}(p_{h}), where P_{*} is defined by
$\left({P}_{}p\right)\left({x}_{1},{x}_{2}\right)={\mathrm{esssup}}_{{x}_{3}\left(,\right)}p\left({x}_{1},{x}_{2},{x}_{3}\right)$
In our numerical examples, we are using the projection P_{*}.
2.4. Numerical Approximation Using Finite Elements
Both elliptic boundary value problems, (2) and (4), are solved with standard finite element approaches [22]: We cover the computational domain Ω (that is, Ω=Ω for the electric potential or $=\stackrel{}{}$ for the temperature) with a uniform Cartesian grid with grid parameter h>0. Let V_{h} be the space of all continuous functions on $\stackrel{}{}$ that are piecewise trilinear (bilinear for (4)). Let Ω_{±,h}, Ω_{pr,h}, and Ω_{ves,h} denote the union of all grid cells that have non-empty intersection with Ω_{±}, Ω_{pr}, or Ω_{ves}, respectively, and let Г_{out,h}, Г_{ves,h}, Г_{±,h}, and Г_{iso,h} be defined as the respective boundary sections. Further, for any topological space X and any set F of functions X→R, any subset $MX$, and any number $cR$ , we use the notation
${F}^{M=c}=\left(fF:F{|}_{M}c\right)$
i. e. F^{M=c} denotes the set of functions in F that take the constant value c everywhere on M.
With this notation, our numerical solution φ_{h} of the potential equation (2) is given by
$\underset{{}_{\mathrm{pr},h}}{\underset{}{}}{}_{h}\left(x\right){v}_{h}\left(x\right)\mathrm{dx}=\underset{{}_{\mathrm{out},h}}{}\frac{n\left(x\right)\left(sx\right)}{{\left(sx\right)}^{2}}{}_{h}\left(x\right){v}_{h}\left(x\right)\mathrm{dA}\left(x\right)$
for all ${v}_{h}{V}_{h}^{\stackrel{}{},h=0}$ and the additional condition
${}_{h}{V}_{h}^{\stackrel{}{}+,h=1,\stackrel{}{},h=1}$
This results in a system of linear equations for the nodal values of φ_{h} that offers a unique solution. The temperature equation (4) is discretized analogously.
2.5. Experimental Setup
2.5.1. Computational Experiments
We performed the patient independent precomputation as described in this paper for a virtual RF applicator that was modelled according to the properties of a CelonProSurge 150-T30™ applicator (Celon AG, Teltow, Germany; cf. Sect. 2.5.2 below). Further, we assumed an electric generator with an inner resistance of 80 Ω, set up to a power of 30 W. In the two-dimensional computation of the temperature distribution, we used a computational domain $\stackrel{}{}$ with an extent of 60×60 mm, covered with a grid of width 0.1 mm. For the three-dimensional computation of the electric potential, we used a domain Ω of extent 38×38×50 mm with a grid width of 0.1 mm in x_{1} and x_{2} directions, but a width of 1 mm in x_{3} direction. The tissue was assumed to denaturate at a temperature of T_{crit} = 50 °C. The material parameters σ, λ, and ν were chosen as follows [14]:
$=0.21\frac{S}{m}\text{}=0.437\frac{W}{\mathrm{Km}}$
and
$v={v}_{\text{0}}{}_{\mathrm{blood}}{c}_{\mathrm{blood}}\mathrm{where}$
${v}_{0}\text{}=0.01765{s}^{1}\text{}\left(\text{relative blood flow rate}\right)$
${}_{\mathrm{blood}}\text{}=1059\mathrm{kg}/{\text{m}}^{3}\text{}\left(\text{blood density}\right)$
${c}_{\text{blood}}\text{}=3850J/gK\text{}\left(\text{blood heat capacity}\right)$
2.5.2. Ex-Vivo Experiments
A first clinical ex-vivo evaluation of our model has been performed. The clinically relevant part of the evaluation is published separately [27]; we focus here on the aspects that are essential in connection with the current publication.
Pieces of fresh porcine livers were ablated using a CelonProSurge 150-T30™ applicator (bipolar, active zone length 30 mm, isolator length 3 mm, radius 0.9 mm). The power unit was set up to p_{setup} = 30W, but note that due to the varying tissue impedance, the effective power was generally below that value. The ablation was performed until a total energy amount of 15 kJ had been induced (where the required duration varied, see Sect. 3.2 below for details). A glass tube, arranged parallel to the applicator and connected to a pump with an adjustable flow rate of saline, has been used to simulate a blood vessel. The experiments were performed at room temperature (25˚C). A number of experiments has been performed that varied in
- the distance d_{ap} of the glass tube from the applicator (5 mm or 10 mm),
- the radius r_{ves} of the glass tube (outer radius between 2.5 mm and 5 mm), and
- the flow rate in the glass tube (between 8.3 ^{cm}/_{s} and 330 ^{cm}/_{s}, which includes the range that has been measured in hepatic blood vessels [15, p. 79], [28, p. 62]).
After each experiment, the liver was cut perpendicular at the center of the applicator. A photograph was taken and imported into a software tool. The lesion border, applicator center, and glass tube location were marked manually with that software tool.
3. RESULTS
3.1. Computational Results
The computations described in Sect. 2.5.1 generated a sparse look-up table containing 526 nodes in the (r_{ves},d_{ap}) plane, where 65 values in the direction α_{nec} were stored. The computation required less than $2\frac{1}{2}$ hours on a four-years-old standard desktop computer, where about 20 minutes of this time were required to compute the electric potential φ and power density p (cf. Sect. 2.3).
Fig. (7) shows a reconstruction of the necrosis area from the data stored in the look-up table for various values of r_{ves} and d_{ap}. We did not consider vessels with a radius below 1 mm (diameter below 2 mm) since these will usually be destroyed by blood clots [29]. As can be seen in the figure, encasement can happen (although only rarely), which justifies the additional consideration of a tube of unaffected tissue around the vessel. On the other hand, in cases where encasement happens, the tube is very thin (fractions of a millimeter). This is confirmed by Fig. (8), which shows a plot of the tube thicknesses (values of Ψ) as functions of d_{ap} for various values of r_{ves}: The tube thickness remains below a clinically significant value for settings in which encasement is imaginable to occur; note that encasement can in particular not occur for ${d}_{\mathrm{ap}}\underset{}{>}9\mathrm{mm}$ since the maximally obtained coagulation radius is about 10 mm.
Since considering this tube does not slow down the computation essentially (neither when creating the tables nor when applying the result to a patient specific data set), we currently decide to retain it.
An interactive software tool, based on the MeVisLab platform (www.mevislab.de), has been developed that is meant to assist the physician in the planning of the intervention. It allows for placement and rearrangement of one or several RF applicators and uses colors and arrows to point out underablated tumor regions. See Fig. (9) for an example obtained from a clinical real CT data set, from which the tumors and vessel trees were segmented and preprocessed with methods presented earlier [19, 20].
In the visualization, every point ξ on the tumor surface has to be checked for its containedness in C, which involves a double loop over the points of the vessel skeleton and the electrode axis. We speed up the computation by instead of the outer loop, computing the projection x^{*} = p_{[a*, b*]}(ξ) of ξ onto the electrode axis and replacing the check for $C$ with the check
$\underset{i=1}{\overset{N}{}}\underset{y\left({a}_{i},{b}_{i}\right)}{}\left({C}_{y,{r}_{i},{r}^{}}^{\left(1\right)}\left({x}^{}\right){C}_{y,{r}_{i},{r}^{}}^{\left(2\right)}\left({x}^{}\right)\right)$
for this x^{*}, which, although it potentially underestimates the necrosis, did apparently not influence the result in the examined cases.
3.2. Ex-Vivo Results
The mean duration of the experiments described in Sect. 2.5.2 was 985 s with a standard deviation of 141 s (n=361). Three example results are shown in Figs. (10-12); note again that the clinically relevant part of the ex-vivo study is published separately [27].
It turned out that between intervention and measurement, a deformation of the tissue and, in particular, a shrinkage of the necrosis took place, which made a quantitative comparison challanging (see figures).
Nevertheless, the evaluation can be used to check the computation results for plausibility. As can be seen from the images, the shape of the necrosis area apparently coincides with the computed results (cf. framed plots in Fig. 7), and so do the total size and, in particular, the distance between the vessel and the necrosis (θ_{tube}).
Additionally, the results give rise to the conclusion that the blood flow rate has essentially no influence on the ablation result [27]. This implies that the blood is not heated up considerably, i. e. the Dirichlet boundary condition for the temperature (4c) is justified. Other authors consider the so-called thermal equilibration length, that is the length of a blood vessel within which the blood is heated up such that the temperature difference to the surrounding tissue is reduced by a factor of exp(–1). According to Creeze and Lagendijk [30], the thermal equilibration length for vessels of the size we are concerned about is essentially above the length scale which the heat of an RF applicator affects, which confirms our observation.
Another point is that the ex-vivo results do not show a significant dependence on the vessel radius [27]. The computational results clearly indicate such a dependence. However, this is not a contradiction, since this dependence mainly appears if r_{ves} and / or d_{ap} are chosen outside the range covered by the ex-vivo experiments (that is, r_{ves} < 2.5mm or d_{ap} < 5mm), see again Fig. (7).
Additional ex-vivo experiments were performed to verify that the glass did not influence the necrosis: We compared experiments with disabled pump to those without any glass tube, not obtaining visual differences [27].
4. DISCUSSION
We presented a concept for a fast estimation of the cooling effects induced by the vascular systems in hepatic RF ablation. The idea is based on patient-independent precomputations, the results of which are stored in a look-up table. Currently, the look-up table establishes a relationship between the radius of a vessel, the distance of an RF applicator from the vessel, and the cooling effect due to this vessel. The cooling effect is expressed as a direction dependent size of the necrosis, together with the distance from the vessel in which tissue remains unaffected.
The results of our computation can be visualized using a tool that allows to place and move applicators interactively and is meant to assist the physician in the intervention planning step. First ex-vivo experiments confirmed some of the computational results.
We would like to emphasize that our visualizations are not meant to replace a thorough intervention planning, nor are they meant to claim that in certain situations RF ablation cannot be successful. It is of course up to the physician’s choice to perform additional steps like a Pringle manoeuvre or chemoembolization to avoid the cooling effect.
In our future work, we plan a number of extensions to our model, including the following:
- We will replace the currently used, partially two-dimensional forward simulation model (Sect. 2.3) with a fully three-dimensional, time-dependent model that also takes state dependent material paramters into account.
- We will examine whether the orientation of the applicator (in relation to the blood vessel) has an essential influence on the result and include this parameter into the look-up table, if necessary. This will introduce additional dimensions in the look-up table, hence requiring essentially more precomputations.
Other subjects of future work are a thorough evaluation, including a verification of the assumption of additive cooling (cf. Sect. 2.1), as well as the creation of look-up tables for a large number of clinically employed applicator and generator types and settings. The validity of our modelling of the blood vessels (constant temperature, same electric conductivity as in the parenchyma) will also be examined.
In view of a more realistic evaluation, it would be useful to employ a setup similar to that described by Lubienski et al. [31].
As has been published earlier [32, 33], the cooling effect of different hepatic vascular systems (that is in particular, the portal venous (PV) versus the hepatic venous (HV) tree) result in considerably different necrosis shapes, see Fig. (13). However, to our knowledge, no forward simulation model for the incorporation of this effect has been proposed by now. Once such a model has been developed, we will incorporate it into our concept. Note that this is straight forward by just preparing independent look-up tables for these two systems (roughly doubling the number of precomputations).
Finally, we would like to remark that our proposed method strongly relies on a precise vessel and tumor segmentation result. The radius information supplied by the currently employed segmentation algorithm [19, 20] is, however, not accurate enough for this purpose [34]. Note that even although the vessel radius is less decisive for the cooling effect (cf. Sect. 3.2), an incorrect vessel radius will (supposed that the vessel centerline has been determined correctly) also essentially distort the information about the distance of both applicator and tumor from the vessel. Once better segmentation procedures have been developed, they can directly be used by our model—recomputation of the look-up tables will not be necessary.
ACKNOWLEDGEMENTS
Part of this research was supported by the German Research Foundation (project Pe199/15-2). We would like to thank Thomas Stein from Celon AG, Teltow, Germany, for supplying precise information about the utilized applicators, as well as Susanne Zentis, Claudia Hilck, and Michaela Jesse from MeVis Medical Solutions for performing the analysis of the CT scans from the example case. We also thank all our colleagues at Fraunhofer MEVIS for hints and helpful discussions and in particular Andreas Weihusen and Christian Rieder for their work on software toolkits for the clinical evaluation.
REFERENCES
[1] | World Health Organization, "“Fact sheet no. 297 ‘cancer”, Feb. 2009. [Online]", Available: http://www.who.int/mediacentre/factsheets/fs297.Accessed: Oct. 5, 2009 |
[2] | P.L. Pereira, "“Actual role of radiofrequency ablation of liver metastases,”", Eur. Radiol, vol. 17, pp. 2062-2070, 2007. |
[3] | S.N. Goldberg, C.J. Grassi, J.F. Cardella, J.W. Charboneau, G.D. Dodd III, D.E. Dupuy, D Gervais, A.R. Gillams, R.A. Kanna, F.T. Lee, T Livraghi, J McGahan, D.A. Phillips, H Rhim, and S.G. Silverman, "“Image-guided tumor ablation: Standardization of terminology and reporting criteria,”", Radiology, vol. 235, pp. 28-739, 2005. |
[4] | D.S.K. Lu, S.S. Raman, D.J. Vodopich, M Wang, J Sayre, and C Lassman, "“Effect of vessel size and creation of hepatic radiofrequency lesions in pigs: assessment of the “heat sink” effect,”", Am. J. Roentgenol, vol. 178, pp. 47-51, 2002. |
[5] | S Lindquist, "“The heat-shock response,”", Ann. Rev. Biochem, vol. 55, pp. 1151-1191, 1986. |
[6] | G Schueller, J Kettenbach, R Sedivy, A Stift, J Friedl, M Gnant, and J Lammer, "“Heat shock protein expression induced by percutaneous radiofrequency ablation of hepatocellular carcinoma in vivo,”", Int. J. Oncol, vol. 24, pp. 609-613, 2004. |
[7] | D.D. Mosser, A.W. Caron, L Bourget, C Denis-Larose, and B Massie, "“Role of the human heat shock protein hsp70 in protection against stress-induces apoptosis,”", Mol. Cell Biol, vol. 17, pp. 5317-5327, 1997. |
[8] | L Frich, K Bjørnland, S Pettersen, O.P.F. Clausen, and I.P. Gladhaug, "“Increased activity of matrix metalloproteinase 2 and 9 after hepatic radiofrequency ablation,”", J. Surg. Res, vol. 135, pp. 297-304, 2006. |
[9] | R Wagenaar-Miller, L Gorden, and L.M. Mastrisian, "“Matrix metalloproteinase in colorectal cancer: Is it worth talking about?,”", Canc. Metastat. Rev, vol. 23, pp. 119-135, 2004. |
[10] | S Mulier, Y Ni, J Jamart, L Michel, G Marchal, and T Ruers, "“Radiofrequency ablation versus resection for resectable colorectal liver metastases: time for a randomized trial?,”", Ann. Surg. Oncol, vol. 15, pp. 144-157, 2007. |
[11] | C Villard, L Soler, and A Gangi, "“Radiofrequency ablation of hepatic tumors; simulation, planning, and contribution of virtual reality and haptics,”", Comput. Meth. Biomed. Eng, vol. 8, pp. 215-227, 2005. |
[12] | T Stein, Untersuchungen zur Dosimetrie der hochfrequenzstrominduzierten interstitiellen Thermotherapie in bipolarer Technik, Fortschritte in der Lasermedizin, vol. 22. 2000. |
[13] | S Tungjikusolum, S.T. Staelin, D Haemmerich, J.-Z. Tsai, H Cao, J.G. Webster, F.T. Lee, D.M. Mahvi, and V.R. Vorperian, "“Three-dimensional finite-element analyses for radio-frequency hepatic tumor ablation,”", IEEE Trans. Biomed. Eng, vol. 49, pp. 3-9, 2002. |
[14] | R Larsen, M Nielsen, and L Sporring, “9th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI)”, Lecture Note in Computer Science. vol. 4191, pp. 380-388, 2006. |
[15] | C Welp, "“Simulationsbasierte Untersuchung der hochfrequenten Thermoablation bei Lebertumoren,”", Ph.D. thesis, Ruhr-Universität Bochum: Bochum, Germany, 2005. |
[16] | J Mooibroek, and J.J. Lagendijk, "“A fast and simple algorithm for the calculation of convective heat transfer by large vessels in three-dimensional inhomogeneous tissues,”", IEEE Trans. Biomed. Eng, vol. 38, pp. 490-501, 1991. |
[17] | T.W. Sheu, C Chou, S Tsai, and P Liang, "“Three-dimensional analysis for radio-frequency ablation of liver tumor with blood perfusion effect,”", Comput. Meth. Biomech. Biomed. Eng, vol. 8, pp. 229-240, 2005. |
[18] | E.J. Berjano, "“Theoretical modelling for radiofrequency ablation: state-of-the-art and challenges for the future,”", Biomed. Eng. Online, vol. 5, no. 24, April 2006.Accessed Oct. 5, 2009 [Online]. Available: http://www.biomedical-engineering-online.com/content/5/1/24 |
[19] | H Bourquain, A Schenk, F Link, B Preim, G Prause, and H.O. Peitgen, "“Hepavision2: a software assistant for preoperative planning in living-related liver transplantation and oncologic liver surgery,”", Proceedings of the 16th International Congress and Exihibition on Computer Assisted Radiology and Surgery (CARS), pp. 341-346. |
[20] | D Selle, B Preim, A Schenk, and H.O. Peitgen, "“Analysis of vasculature for liver surgical planning.,”", IEEE Trans. Med. Imaging, vol. 21, pp. 1344-1357, 2002. |
[21] | H Bungartz, and D Dirnstorfer, "“Multivariate quadrature on adaptive sparse grids,”", Computing, vol. 71, pp. 89-114, 2003. |
[22] | C Johnson, “Numerical Solution of Partial Differential Equations by the Finite Element Method. Dover Publications: UK, 2009. |
[23] | S.N. Antontsev, and M Chipot, "“The thermistor problem: Existence, smoothness, uniqueness, blowup,”", SIAM J. Math. Anal, vol. 25, pp. 1128-1156, 1994. |
[24] | H.H. Pennes, "“Analysis of tissue and arterial blood temperatures in a resting forearm,”", J. Appl. Physiol, vol. 1, pp. 93-122, 1948. |
[25] | P Deuflhard, M Weiser, and M Seebaß, "“A new nonlinear elliptic multilevel fem applied to regional hyperthermia,”", Comput. Vis. Sci, vol. 3, pp. 115-120, 2000. |
[26] | S Arrhenius, "“Über die Reaktionsgeschwindigkeit bei der Inversion von Rohrzucker durch Säuren,”", Z. Phys. Chem, vol. 4, pp. 226-248, 1889. |
[27] | K.S. Lehmann, J.P. Ritz, S Valdeig, V Knappe, A Schenk, A Weihusen, C Rieder, C Holmer, U Zurbuchen, P Hoffman, H.O. Peitgen, H.J. Buhr, and B.B. Frericks, "“Ex situ quantification of the cooling effect of liver vessels on radiofrequency ablation,”", Langenbeck’s Archiv. Surg, vol. 394, pp. 475-481, 2009. |
[28] | W Wermke, "Sonographische Differenzialdiagnose Leberkrankheiten", Deutscher Ärzte Verlag GmbH, 2006. |
[29] | T De Baere, A Denys, B.J. Wood, N Lassau, M Kardache, V Vilgrain, Y Menu, and A Roche, "“Radiofrequency liver ablation: Experimental comparative study of water-cooled versus expandable systems,”", Am. J. Roentgenol, vol. 176, pp. 187-192, 2001. |
[30] | J Crezee, and J Lagendijk, "“Temperature uniformity during hyperthermia: the impact of large vessels,”", Phys. Med. Biol, vol. 37, pp. 1321-1337, 1992. |
[31] | A Lubienski, R.G. Bitsch, K Lubienski, G Kauffmann, and M Duex, "“Radiofrequency ablation (rfa): development of a flow model for bovine livers for extensive bench testing,”", Cardiovasc. Interv. Radiol, vol. 29, pp. 1068-1072, 2006. |
[32] | C Bangard, A Gossmann, H.U. Kasper, M Hellmich, J Fischer, A Hoelscher, K Lacker, and D.L. Stippel, "“Experimental radiofrequency ablation near the portal and the hepatic veins in pigs: differences in efficacy of a monopolar ablation system,”", J. Surg. Res, vol. 135, pp. 113-119, 2006. |
[33] | B.B. Frericks, J.P. Ritz, T Albrecht, S Valdeig, A Schenk, K.J. Wolf, and K Lehmann, "“Influence of intrahepatic vessels on volume and shape of percutaneous thermal ablation zones: In vivo evaluation in a porcine model,”", Invest. Radiol, vol. 43, pp. 211-218, 2008. |
[34] | J Drexl, V Knappe, H.K. Hahn, K.S. Lehmann, B Frericks, H Shin, and H.O. Peitgen, "“Accuracy analysis of vessel segmentation for a litt dosimetry planning system,”", In: T.M. Buzug, and T.C. Lueth, Eds., Perspective in Image-Guided Surgery. World Scientific, 2004. |