Fast Estimation of the Vascular Cooling in RFA Based on Numerical Simulation.

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.

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.

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.

Quantitative Description of Necrosis

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 rves of the vessel and on the distance dap 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 rnec 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 rves, the distance dap of the applicator from the vessel, and the resulting necrosis size rnec, 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

see also Fig. (). For a precise definition we refer the reader to Sect. 2.3.

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

see Fig. (), where we approximate the tube by a circular cylinder.

Full Model

Now, let the complete vascular system be given as a set of straight vessel segments with end points a, b ∈ R3 and radii r > 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. ). 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. (). 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 ∈ R3 to b ∈ R3 and by ∠abc the angle between the lines and . Additionally, by d (a,b) we denote the distance of a ball of radius r around a from a ball of radius s around b, i. e.

Let ξ ∈ R3 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 ,b] 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

We therefore define

so that the set 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

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. ().

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 rves, dap and αnec , the values Ф(rves, dap, αnec) and Ψ(rves, dap) can be computed, but this computation is too slow for clinical application. The idea is to approximate Ф and Ψ by some functions and , respectively, with the following properties:

The construction of and (for a given tolerance) requires preferably few evaluations of Ф and Ψ.

Once and 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 and 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 . We emphasize that no special properties of any particular forward simulation model are utilized here.

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 Ω ⊂ R3 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 x3 direction and have an x2 coordinate of 0, cf. Fig. (). The center of the applicator’s active zone is placed in the origin; the x1 position of the blood vessel is chosen to be negative and such that the distance between blood vessel and applicator matches the prescribed value dap. 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

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]

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 where H1 = H1,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]

and the remaining quantities are properties and settings of the electric generator, namely its inner resistance RI and the setup power psetup. (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 , and if we extend p to Ω by setting p|Ω =0, we have

Performing these computations in two (rather than in three) space dimensions would result in a mismatch of units: ptotal 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 Rtis cause to have a wrong unit as well, and finally the equation for peff 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 denote the intersection of Ω and the x1 - x2 plane as a subset of R2, i. e. (see again Fig. ) and analogously for Ωves, Ωpr, and Гout. Further, we let be given by 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]

where λ>0 is the thermal conductivity, Tbody 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 .

Once T is known, we let Ωcoag ⊂ Ω denote the set of coagulated tissue,

(where Tcrit is a critical temperature above which the tissue is assumed to be destroyed).

Finally, we define the necrosis size rnec 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 H1 is in general not defined, the set Ωcoag is only defined up to sets of measure zero, so that the precise definition of rnec and θtube is given by

where e = (cos α, sin α)and xves is the position of the vessel’s centerline in Ω (i. e. xves = (–rves–dap,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

where 2ℓ is the length of the applicator’s active zone and is a non-negative function with compact support in [-ℓ,ℓ] and the property that . It is straight forward to prove that P is well-defined, but the norm does not remain bounded for ε→0. In other words,

does apparently not provide a mapping , 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 of the power density p (see Sect. 2.4 below) is an element of L∞(Ω), so that P0(ph) is well-defined, and so is even P*(p), where P* is defined by

In our numerical examples, we are using the projection P*.

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 for the temperature) with a uniform Cartesian grid with grid parameter h>0. Let Vh be the space of all continuous functions on that are piecewise trilinear (bilinear for (4)). Let Ω±,, Ωpr,, and Ωves, denote the union of all grid cells that have non-empty intersection with Ω±, Ωpr, or Ωves, respectively, and let Гout,, Гves,, Г±,, and Гiso, be defined as the respective boundary sections. Further, for any topological space X and any set F of functions X→R, any subset , and any number , we use the notation

i. e. F 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

for all and the additional condition

This results in a system of linear equations for the nodal values of φ that offers a unique solution. The temperature equation (4) is discretized analogously.

Experimental Setup

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 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 x1 and x2 directions, but a width of 1 mm in x3 direction. The tissue was assumed to denaturate at a temperature of Tcrit = 50 °C. The material parameters σ, λ, and ν were chosen as follows [14]:

and

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 psetup = 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 dap of the glass tube from the applicator (5 mm or 10 mm),

the radius rves 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.

RESULTS

Computational Results

The computations described in Sect. 2.5.1 generated a sparse look-up table containing 526 nodes in the (rves,dap) plane, where 65 values in the direction αnec were stored. The computation required less than 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. () shows a reconstruction of the necrosis area from the data stored in the look-up table for various values of rves and dap. 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. (), which shows a plot of the tube thicknesses (values of Ψ) as functions of dap for various values of rves: 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 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. () 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(ξ) of ξ onto the electrode axis and replacing the check for with the check

for this x, which, although it potentially underestimates the necrosis, did apparently not influence the result in the examined cases.

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. (); 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. ), 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 rves and / or dap are chosen outside the range covered by the ex-vivo experiments (that is, rves < 2.5mm or dap < 5mm), see again Fig. ().

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].

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. (). 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.