## Abstract

The treatment of cancerous tumours in the liver remains clinically challenging, despite the wide range of treatment possibilities, including radio-frequency ablation (RFA), high-intensity focused ultrasound and resection, which are currently available. Each has its own advantages and disadvantages. For non- or minimally invasive modalities, such as RFA, considered here, it is difficult to monitor the treatment *in vivo*. This is particularly problematic in the liver, where large blood vessels act as heat sinks, dissipating delivered heat and shrinking the size of the lesion (the volume damaged by the heat treatment) locally; considerable experience is needed on the part of the clinician to optimize the heat treatment to prevent recurrence. In this paper, we outline our work towards developing a simulation tool kit that could be used both to optimize treatment protocols in advance and to train the less-experienced clinicians for RFA treatment of liver tumours. This tool is based on a comprehensive mathematical model of bio-heat transfer and cell death. We show how simulations of ablations in two pigs, based on individualized imaging data, compare directly with experimentally measured lesion sizes and discuss the likely sources of error and routes towards clinical implementation. This is the first time that such a ‘loop’ of mathematical modelling and experimental validation *in vivo* has been performed in this context, and such validation enables us to make quantitative estimates of error.

## 1. Introduction

The minimally invasive nature of thermal ablation (the heating of tissue to induce death) makes it a common choice for the treatment of primary hepatocellular cancerous tumours in the liver when the prognosis of surgical resection techniques is poor. Unfortunately, the success rate drops off rapidly with decreased experience: for liver tumours, clinicians with 2–4 years experience achieve 89 per cent survival rates at 2 years, significantly above the 46 per cent survival rate managed by those with fewer than 2 years of experience [1]. This is largely owing to the complex relationship between thermal power and final lesion size, which is very strongly affected by localized perfusion and the presence of nearby large cooling vessels, which act as significant heat sinks in the liver. It is thus very difficult to tailor the heating protocol to the individual patient.

During thermal ablation, such as radio-frequency ablation (RFA), which uses a high-frequency alternating current to induce heating, visualization of the progress of the resulting thermal coagulation is often unavailable. As a result, mathematical modelling plays a vital role in the determination of *in vivo* temperatures to plan and optimize treatment protocols [2,3]. It would be extremely valuable to have a computational tool that, based on individual patient imaging data, could be used to optimize the heating protocol in advance of the treatment. This could also be used as a training tool for the less-experienced clinicians to improve outcome rates where they are at their poorest level.

However, the development of such a tool is extremely difficult because an extensive validation is required before it can be used with confidence in a clinical setting. Such validation requires predictions to be made of the size of the lesion, based solely on routinely available imaging data for either human or animal models, and then compared with the status of the liver tissue after treatment. This requires a combination of mathematical modelling, computational efficiency, accurate image analysis and carefully designed experimental studies; the need for such a team to be gathered to ‘close the loop’ between prediction and reality is a significant obstacle to achieving this.

In this paper, we present our progress towards this and show the results that we have obtained on animal models. This has been performed through the development and implementation of a multi-scale mathematical model that can be used to simulate the ablation process, based on patient-specific data and images. Using the image analysis, we are able to perform a quantitative validation of the model predictions of lesion size compared with those acquired through histological examination. These are presented and discussed alongside plans to extend the work into human models and hence into clinical practice. This would be the first such simulation tool kit for this therapy in the liver, and has the potential to make a substantive impact on clinical outcomes, as well as providing valuable pointers for future successful work in other organs.

## 2. Mathematical model

A schematic of the mathematical model is shown in figure 1, illustrating how the different parts of such a model connect. The two major components are the bio-heat model, which uses the heating and blood flow field to predict the tissue temperature (which varies both temporally and spatially), and the tissue state model, which requires the tissue temperature profile to predict the state of the tissue (alive or dead). This status in turn affects the tissue parameters, such as thermal conductivity (which may also be temperature dependent), and hence influences the bio-heat model. There is thus considerable interaction between the different sections of the model, which exhibits many nonlinearities.

In order to predict the outcome of a heating protocol, four things are thus required: the heating profile, the blood flow field, a bio-heat model (with parameter values) and a tissue state model (with parameter values). For a review of the literature, see Payne *et al*. [3], where the practical difficulties of choosing models and parameter values are discussed in more detail. In this paper, we describe the practical implementation which we have adopted for our particular application, that of RFA in the liver.

### (a) Bio-heat model

The classical bio-heat equation, the Pennes model, has been widely used in the literature to model the electrical–thermal heating process during the ablation procedure [4–6]. However, the convective heat-transfer term between tissue and blood in the equation is oversimplified, assuming the blood to be a volumetric heat sink and one that is uniformly distributed throughout the tissue [7]. Owing to the absence of an equation to model variations in blood temperature, large deviations have been found between Pennes model predictions and the absolute temperature field in vascularized tissue [8,9]. In order to overcome these shortcomings, many subsequent modifications have been proposed, including those by Wulff [10] and Klinger [11], who considered the local blood mass flux in order to account for the blood flow direction, but continuous equilibrium between tissue and blood temperatures is assumed. Again, the validity of this thermal equilibrium assumption has been questioned under conditions of thermal ablation [12]. Inaccurate estimation of blood cooling effects can greatly affect the prediction of tissue temperature and thus the necrosis zone after the intervention, especially for highly blood-perfused organs such as the liver.

Therefore, in this project, we have proposed a new split-volume bio-heat equation to model the heat transfer during RFA of liver tissue [13]. We adopt the same fundamental physical approach as that of Chen & Holmes [14], who considered each elemental volume to comprise a tissue fraction and a blood fraction. Separate equations are postulated based on conservation of energy for the tissue and blood phases, respectively,
2.1and
2.2where *ρ*, *c*, *k*, *H*, *ε*, **v**, *σ*, *E* and *T* denote density, specific heat capacity, thermal conductivity, local volumetric heat-transfer coefficient, blood volume fraction, blood flow velocity field, electrical conductivity, electric field and temperature, respectively, and subscripts ‘t’ and ‘b’ refer to the tissue and blood phases, respectively. The tissue and blood temperatures are the dependent variables to solve, for which the blood flow velocity field and electric field must be prescribed (as described below). In these equations, the first term on the right-hand side is the radio frequency heating source term, the second the heat-conduction term and the third the tissue/blood convective heat exchange term. The bulk flow term appears on the left-hand side of the blood bio-heat equation as part of the total temporal derivative. The heat exchange term between the blood and the tissue represents convective heat transfer via the vascular wall, which can be obtained through Newton's cooling law,
2.3where *s*_{b} and *h*_{b} are the volume-averaged specific surface area and convective heat-transfer coefficient, respectively. For a bundle of vascular tubes of radius *R*, we can estimate that *s*_{b}=2*V* _{b}/*RV* =2*ε*/*R* and that *h*_{b}=*k*_{b}*Nu*/2*R*, such that *H*=*εk*_{b}*Nu*/*R*^{2} [15], where *Nu* is the Nusselt number, set to be 4.93 for these conditions [16], and *V* denotes volume. This illustrates clearly that the local volumetric heat-transfer coefficient depends on both the blood volume fraction and the vessel radius.

Differing from traditional bio-heat models, the split-volume bio-heat model incorporates the local thermal non-equilibrium between the blood and the peripheral tissue and uses two separate equations to represent the thermodynamics of the blood and solid tissue phases. When used to examine the separate contribution of different parts of the vasculature, the model behaviour is shown to be dependent on a non-dimensional number, *γ*, termed the thermal significance coefficient, which is a measure of the relative heat flow between the blood and the tissue volume phases, as described completely in Peng *et al*. [13]. For vessels in the highest generation of the vasculature such as the big arteries and veins, *γ* is found to be far larger than unity, indicating that the blood in these large vessels essentially holds a constant temperature, in which case, the split-volume bio-heat model can be simplified into the Pennes model; on the other hand, for small vessels in the bottom generations of the vasculature such as arterioles, capillaries and venules, *γ* is much smaller than unity, suggesting that the blood in small vessels is in continuous equilibrium with tissue temperature, in which case, the model is equivalent to the Klinger model. However, most of the temperature equilibration occurs as the blood travels via the middle generations of the vasculature such as terminal artery branches, under which circumstance *γ* is of order unity and hence the split-volume model cannot be further simplified.

In this project, large vessels (diameter greater than 10 mm) can be individually reconstructed from a vessel-enhanced X-ray computed tomography (CT) image, as outlined below, and thus their cooling effects are determined by calculating the exact convective heat transfer on their vessel wall [6]. The cooling effects of very small vessels (diameter less than 0.1 mm) can be considered to be of negligible importance, on the basis that their contribution to vasculature cooling is smallest [13]. Vessels in the middle generations of the vasculature, such as the terminal artery branches, are then modelled as a porous media flow field, as outlined below, and their cooling effects estimated by the split-volume thermal model.

### (b) Cell death model

Once the tissue temperature has been modelled, the second major component of the mathematical model is that which governs the relationship between tissue temperature and tissue state. For this, we have proposed a novel three-state model, outlined in more detail by O'Neill *et al*. [17], extending the existing models that are predominantly based on two states (alive and dead). This three-state model includes an intermediate compartment between the fully alive (*A*) and dead (*D*) compartments, which we term the vulnerable (*V*) state,
2.4where transitions between the states are governed by (temperature-dependent) forward and (temperature-independent) backward rate constants, and *A*, *V* and *D* represent the fraction of cells within each state at each point in the tissue. While most of the literature (see [18] for a very thorough review) uses Arrhenius-based rate coefficients, as proposed by early workers in this area, such as Henriques and Moritz for cutaneous skin burns [19–21], the nature of this exponential form, where the transition rate is proportional to the exponential of the activation energy divided by temperature, results in a highly sensitive constant of proportionality where reported values range over tens of orders of magnitude [18].

Our model is similar to the two- and three-state models proposed by Breen *et al*. [22] and Szasz & Vincze [23], and it has been shown by O'Neill *et al*. [24] that with certain mathematical manipulations, it can be related to the physiologically based model with infinite states as proposed by Jung [25–27]. The biological implication of this is that the vulnerable compartment accounts for cells that are at elevated risk of critical damage.

For simplicity, all complex stages of the path to cell death are incorporated into a single-damage mechanism; thus, as the processes from alive to vulnerable and from vulnerable to dead are taken to be two stages of the same result, a single forward rate constant, *k*_{f}, is used. This rate constant is nonlinear, being proportional to the combined proportion of cells in the vulnerable and dead compartments (1−*A*),
2.5where is a scaling constant and *T*_{k} is a parameter that sets the rate of the exponential increase with temperature. The backward rate constant, simulating any self-healing, is set to a value that is invariant with temperature, as this was found to provide a good fit to the experimental data.

Tailor-made heating experiments for model parameterization were carried out at constant temperatures on cultures made up of two cell lines: human liver hepatocellular carcinoma (HepG2) and human lung fibroblast (MRC-5). Two pure cultures and three co-cultures (75/25, 50/50 and 25/75) were cultured overnight at 37^{°}C. For each time and temperature combination, the culture medium was replaced with a preheated medium and then incubated in a preheated heating cabinet. The temperatures were 55^{°}C, 65^{°}C, 75^{°}C, 85^{°}C and 100^{°}C, and the heating times were 300, 600 and 900 s. Post-heating, the medium was replaced with another medium at 37^{°}C and spiked with 20 μl methanethiosulphonate (MTS) reagent. Fluoroscopy measurements were taken at 2, 26 and 50 h after heating, the readings signifying the metabolic activity of the living cells, i.e. the relative proportions of alive and dead cells. Data were normalized to a control value at 37^{°}C and there were five or six repeats of each plate, which had four separated wells. The data were fitted to the model giving average r.m.s. errors between the model predictions and the data of 1.40 per cent. This compares favourably with the s.e. of 3.39 per cent for the experimental data.

A second part of the model for cell death, which accounts for delayed cell death occurring post-heating after the tissue has returned to normothermia, has also been proposed and validated against the experimental data [17]. Crucially, the data show that there appears to be a threshold for cell viability, estimated to be 80 per cent, such that if the local viability is reduced beyond this level, the cells will inevitably progress to a dead state eventually. This is of particular value computationally, as the simulations need not be run past the time at the end of heating.

Of further significant benefit to the multi-scale simulation of RFA treatments is that the mixed co-culture experiments have highlighted different thermal resistances of different cells. We note that this will be of considerable importance in future refinement of the model for patient-specific treatment planning, where conditions may vary widely.

### (c) Electrical heating model

At the frequencies employed in RFA (300 kHz–1 MHz) and within the area of interest (it is known that the electrical power is deposited within a small radius around the active electrode), the tissues can be considered purely resistive as the displacement currents are negligible. Using a quasi-static approach, the distributed heat source *q* (Joule loss) is given by [28]
2.6where **J** is the current density and **E** is the electric field intensity. The values of these two vectors are evaluated using Laplace's equation,
2.7where *V* _{PD} is the voltage. Joule heat can thus be written as
2.8In the literature, the voltage is usually set to be a constant number on the surface of the radio-frequency (RF) electrode probe needle and the potential field generated around the probe is solved numerically, using equation (2.7) (e.g. Panescu *et al*. [5]). This approach, though theoretically accurate, suffers from a high computational cost because of the required fine volume mesh on the surface of the very thin probe tips. Moreover, the exact position of the whole electrode probe is not always available in CT imaging [3], and the calculated Joule heat has been found to be very sensitive to the probe position [29]. Therefore, it is highly desirable that the deposited Joule heat during RFA can be estimated without prior knowledge of the whole electrode geometry in the tissue.

In practice, the electric field strength generated by an RF electrode probe is only sufficiently high to cause significant direct heating very close to the electrodes, thermal conduction contributing most of the tissue heating at distances farther from the electrode [30]. This raises the possibility of enclosing the thermal energy generated locally by the electrode inside a few voxels containing the electrode probe. Thus, in our electro-thermo simulation, we use a few point sources to represent the power delivered by the electrode, based on the estimation from an analytical analysis of charge distribution and electrical field of the radiofrequency interstitial tumor ablation (RITA) probe used in the animal model experiments here (StarBurst Radiofrequency Ablation, AngioDynamics; www.angiodynamics.com).

As the full details of the analytical study will be presented separately, a brief overview is given here. The analysis starts by solving the charge distribution on a line conductor of finite length, a simplified geometry of the electrode probe. Coulomb repulsion pushes charges out towards the ends, just as the charge on a solid conductor flows to the surface. However, no explicit analytical solution of the charge distribution has been found in the literature; in fact, it has been suggested that the problem itself may be ill-posed, the answer depending on the particular model used to represent the conductor [31].

In our approach, based on the previous study of Andrews [32], we assume that the charge density is uniform apart from at the ends. The radio between the charge density at two ends of the conductor and the middle part can be obtained using a moment method of formulation, given as
2.9where *γ*≈0.5772 is Euler's constant, *L* is the conductor length, *r* is the conductor radius and the end of the conductor Δ*l* is assumed to be 1/*N* of the whole length of the conductor. With the charge density ratio, the joule heat generated by the probe end and the middle part is given as (assuming the total deposited power to be *P*)
2.10and
2.11The RITA probe needle used in the project has a diameter of 0.5 mm. The needle length varies from 1 to 4 cm and hence we use the average, *L*=2.5 cm. The selection of section number *N* is, however, somewhat arbitrary but must satisfy a preliminary condition *N*=*L*/Δ*l*≪*L*/*r* [19]. We choose *N*=10 here, thus obtaining a power ratio of *q*_{end}/*q*_{mid}≈2.5.

Thus, as shown in figure 2, 21 electrical sources were placed in the simulation in the following positions: tip ends (9), tip middles (9), tip intersection (1), probe shaft's middle (1) and probe shaft's end (1). The magnitude of power delivered at the tip ends and the probe shaft's end will be much higher than that of the rest of the electrical sources. Assuming the whole power input is *P*, the power distribution will be: tip ends (9×2.5/36*P*), tip middles (9×1/36*P*), tip intersection (1/36*P*), probe shaft's middle (1/36*P*) and probe shaft's end (2.5/36*P*).

Although this approach is a simplification to the actual heating distribution, the improvement in computational time required is a considerable bonus, given the need to perform simulations rapidly in a clinical setting. Examination of the model predictions with this model shows that the heating distribution is ‘smoothed out’ by the conduction of heat within the tissue and convection within the blood stream, thus the particular details of the heating point sources appear to have only a limited effect on the final lesion size. The most important factor appears to be that the heating is distributed along the probe tips and shaft.

### (d) Flow model

The continuum porous media flow field can be constructed as an extension of the existing CT-visible vasculature, by setting porous media flow sources and sinks at boundaries of the vessel tree to mimic the flow of blood from the large vessels to those in the middle generations of the vasculature (modelled using a porous medium approach). This approach incorporates the patient-specific information of individual vascular geometries and hence better represents the vascular-mediated cooling. As a simplified case, this model was used to analyse the cooling effects of a CT-visible vessel and its unregistered branches [33]. Strong directional effects of perfusion on the tissue temperature distribution have been found, which is confirmed by recent experiment results [34]. Therefore, the same model equations are applied to the real vascular geometries in the liver, as presented in this paper.

Using the theory of porous media, Darcy's law states that the flow velocity is the gradient of pressure,
2.12where *κ* is the permeability of the medium and *μ* is the dynamic viscosity of the fluid. The pressure field, *p*, can be evaluated using Laplace's equation based on continuity of blood flow velocity,
2.13To solve equation (2.13), boundary conditions must be specified on the organ boundary and the CT-constructed vessel tree. In our model simulation, mixed boundary conditions are used by setting either flow velocity *u*_{b} or pressure value *p* on all boundary points. On the organ boundary, a Neumann boundary condition is used, given as *u*_{b}=0; while on the CT-constructed vessel tree, a Dirichlet boundary condition is used: a small positive pressure is applied on the portal vein vessel wall and a small negative pressure is applied on the hepatic vein vessel wall. This drives the porous medium flow, which mimics the middle generations of the vasculature, as described above.

## 3. Image analysis

As outlined in the mathematical model in §2.4, we now examine how the geometric models were generated. The data included images from two pigs, both during and after RFA. The images were acquired at arterial, portal venous and hepatic venous phases and with forced breath-hold techniques to minimize breath-related motion. The data were imaged at the Medical University of Graz using a 320-line Aquilion ONE CT scanner (Toshiba Medical Systems, Otawara, Japan). The images were scanned with a matrix size of 512×512×320 and a resolution of 0.4 and 0.5 mm for in-plane and out-of-plane, respectively.

In the image analysis, we extracted the geometric models of the liver, hepatic artery, portal vein, hepatic vein and location of RFA needle tips. In this paper, healthy pigs were used and thus tumour segmentation was not needed. The locations of needle tips were extracted from intra-operative volumes (imaged during the RFA treatment), whereas anatomic structures were extracted from post-operative volumes (imaged after the RFA treatment). With registration, each structure was guaranteed to be spatially aligned in the reference frame.

The proposed image analysis framework has five steps, outlined in detail below: (i) registration of all data to the selected reference frame, (ii) segmentation of liver outline, (iii) segmentation of vessels, (iv) segmentation of the induced lesion, and (v) manual processing and surface meshing.

### (a) Registration

Not all deformations could be removed with the forced breath-hold technique, so an additional registration was used. We used in-house software similar to Rueckert *et al*. [35]. Such approaches have been previously used in liver motion estimation [36]. In the registration, one of the volumes was assigned as a target, and all other volumes denoted as source volumes were registered one by one with the target; first by rigid and then by non-rigid transformation. The non-rigid transformation was parametrized as a cubic B-spline free-form deformation model with control points positioned in a regular lattice [35,36].

Normalized mutual information (NMI) was used to measure the similarity between the volumes [37]. In optimization, the inverse of the NMI between the target and the source volume was minimized with respect to the transformation parameters. We used simplex minimization for the rigid model and conjugate gradient minimization for the non-rigid model [38]. To speed-up the computation, we used hierarchical coarse-to-fine optimization. We used a two-level image pyramid, where the coarse-level resolution was 2.0 mm and the fine resolution was 1.0 mm in each direction. Similarly, the control point distance was hierarchically refined from 40 to 20 mm and finally to 10 mm.

### (b) Liver segmentation

We have developed a new hybrid segmentation method for liver outline. The method is based on landmarking followed by direct and indirect liver segmentation. In landmarking, the goal was to extract the structures that are easy to segment and use them to limit the segmentation in the internal part of the abdominal cavity. Our approach is similar to the landmarking presented by Rangayyan *et al*. [39]. First, the contrast-enhanced CT image was divided in two parts as follows. Air was segmented and components connecting the image borders were marked as a background. The background was morphologically opened to remove external artefacts, the largest internal component was marked as a body, and all other tissues were marked as non-body. Subsequently, the body was divided into thoracic and abdominal cavities based on the location of the diaphragm. We segmented the lung lobes and fitted two quadratic surfaces representing the diaphragm through the bottom of the lobes [39]. The surfaces were combined and the thoracic cavity above the diaphragm was removed.

Finally, we divided the abdomen into internal and peripheral parts. The skeleton was segmented with the region-growing algorithm. We then formed an intact skeletal armour by closing the skeleton in the head–feet direction. The previously extracted non-body tissue and the armour were assigned in the same class and morphologically closed. This removed all the peripheral tissue located inside the body but external to skeletal bones.

In indirect segmentation, the internal part of the abdominal cavity was further processed by iteratively removing non-liver tissue. First, fat and internal air were segmented and morphologically opened with large structuring elements. This process divided the abdominal cavity into separate organs, the largest of which is the liver.

In direct segmentation, the intensity distribution of the healthy liver tissue was estimated. The indirect segmentation above provided an initial liver mask for intensity distribution estimations. The liver was segmented with region growing. Holes, e.g. vessels, were then closed with morphological operations. The process was repeated to all registered post-operative volumes: arterial, portal venous and hepatic venous phase volumes. This results in six different liver segmentations for a single liver. To form a consensus result, a voxel-wise majority of voting was applied to indirect and direct segmentations.

### (c) Vessel segmentation

The hepatic vessels were extracted employing the tried and true skeleton-based approach first introduced by Selle *et al*. [40], further developed with multi-scale and matched filter approaches, intensity ridge-tracking and mathematical morphology.

Matched filters refer to Hessian-based vessel enhancement filters, which characterize local scale-dependent second-order variations around a shape primitive, and thus the likelihood of a vessel, by the eigenvalues of the Hessian matrix of a Gaussian-filtered image. The standard deviation of the Gaussian kernel is chosen according to the desired vessel size, and filter responses with several vessel sizes are merged into one [41]. To make computation practical, multi-scale Hessian convolution was implemented using the image-pyramid approach. The result both enhances vessels and removes noise, so no additional pre-processing is required.

The hepatic vessels were segmented in two phases—the first being wavefront region growing with user initiation. Wavefront implementation supports high-level control over propagation, e.g. topological constrains, and it resolves large- to medium-sized vessels using a single threshold. The smallest vessels were segmented in the second phase, which is a ridge-oriented alternative to earlier locally adaptive schemes, but which differs from the usual merging and reconnection schemes. In this iterative scheme, the initial segmentation is extruded graph-wise using a watershed transform [42]. It can be thought of as propagating both at the watershed line (intensity ridge) and towards the watershed minima, from which the former was extracted using heuristic graph analysis. The method is inherently capable of growing over small gaps and extracting structures down to unit voxel thickness. Finally, a marker-controlled watershed transform for vessel boundaries was performed to complete the segmentation of small vessels.

In post-processing, the segmented vessels were converted into piecewise linear skeletons using a top-down minimum-cost traversal balanced by two distance fields. Conversion to the skeletal domain reveals the vessel tree structure for qualitative analysis, and the transform is approximately invertible, thus it can be used as an intermediate step within the segmentation pipeline. In our previous work [43], this method was proven to be high performing and capable of extracting 97 and 75 per cent of the vessels equal to or above the radius of 2.0 or 1.0 voxel units, respectively. With regard to sensitivity, the results showed improvements over the existing methods.

### (d) Lesion segmentation

The applied lesion segmentation process was similar to direct segmentation of the liver. The intensity distribution of the lesion was estimated and a region-growing algorithm was applied. The result was morphologically opened to remove overflow to adjacent structures. The process was repeated for all registered post-operative volumes. The voxel-wise majority voting was applied for final segmentation.

### (e) Manual refinement and meshing

After the automatic and semi-automatic segmentation, manual refinement of the segmentation results was performed. For liver and lesions, slice-by-slice correction, with ITK-SNAP [44] (www.itksnap.org), was performed. For vessels, we used in-house tools developed for semi-automatic pruning of a vessel tree in the skeletal domain. The needle tips were extracted by hand using intra-operative series. Finally, the geometric models of anatomic structures were produced from image data using the marching cubes approach [45].

### (f) Conclusions

An example of the reconstructed liver geometry is shown in figure 3, illustrating the different components of the vascular tree and the lesion. The high quality of the reconstruction enables the mathematical model to be implemented with confidence, as described in the subsequent section.

## 4. Numerical model implementation

In this section, we describe the practical implementation of the mathematical model outlined in §2, using the physical geometries generated, as outlined in §3. In particular, we focus on the approaches taken to minimize computational expense, as this will be a critical requirement if treatment protocols are to be optimized and hence used in clinical practice.

### (a) Finite-element method discretization of the coupled differential equations

The numerical implementation of the coupled differential equations is performed using the finite-element method, resulting in a spatial discretization of the computational domain. The tetrahedral mesh used for the computation is obtained by a process described at the end of this section.

Each equation is integrated over time using a backward difference formula of order 2, which in the case of the bio-heat equations amounts to expressing the diffusive, convective and forcing terms with respect to the new value of the temperature. This method (implicit in time) has the advantage of being stable and reads, for a second-order discretization of the time derivative,
4.1where the superscript ‘*n*’ denotes the approximation at time *n*Δ*t*, Δ*t* being the time step of the simulation and the case for *n*=0 is treated with a first-order discretization of the time derivative. The full linear system obtained from the discretization of the complete model is solved iteratively with a bi-conjugate gradient method and preconditioned with an incomplete lower and upper triangular matrix factorization.

### (b) Implementation with the ELMER code

The core finite-element procedures are based on the open source multi-physics code Elmer [46]. The existing solvers written in Fortran 90 for basic models (advection–diffusion equation, Poisson equation) have been modified to deal with our particular model. The three resulting solvers (bio-heat equations, porous flow and cell death model) are weakly coupled and called iteratively within each time step.

### (c) Bio-heat equations implementation

The coupled blood and tissue equations are solved by strong coupling through the heat transfer between the two domains. They are nonlinear owing to the linear dependency of thermal conductivity on temperature and are linearized with a fixed-point algorithm (see table 1 for the values of model parameters used here). Assuming that we are dealing with only mild variations of the properties with temperature, the Picard linearization can be used as an iterative method. The heat transfer is through convection and diffusion, therefore the finite-element implementation of the heat equation uses the Galerkin least-squares stabilization to prevent numerical oscillations. In order to model the various heat-transfer conditions at vessels walls, both fixed and convective boundary conditions are implemented for tissue and blood temperatures. A *T*=37^{°}C condition allows a strong cooling effect from blood, whereas the convective transfer produces a relaxation of the preceding condition. The convective transfer can be described by the following relation:
4.2where *α* is the heat-transfer coefficient and *n* the outer normal to the boundary. *T* is measured in degrees Celsius, relative to body temperature (assumed to be constant at 37^{°}C). Equation (4.2) is applied at the walls of the (large) individually registered vessels, as described earlier. The heat-transfer coefficient is expressed as a function of the blood conductivity and the vessel radius. A fixed *T*=37^{°}C condition is applied on the outer wall of the liver to represent the heat-adjacent organ heat sink.

### (d) Joule heating modelling

As the joule heating is proportional to the square of the gradient of the electric potential distribution throughout the liver tissue, we would normally have to solve the Poisson equation with a fine resolution near the needle surface in order to account for the steep gradients. One way to improve the solver efficiency is to approximate the joule heating by lumping the heat generation distribution into specified points on the needle (see §2.3 and figure 2). The sum of the heat generated adjacent to the needle surface can then be added into the heat equation directly as power source terms at these specific points and smoothed with a Gaussian function. The total electric power, which is time-dependent, and tracked by the StarBurstXL generator, is then shared between all the smoothed point sources. In order to approximate the sharp drop in tissue electric conductivity owing to the presence of vapour when the temperature goes above 100^{°}C, we set the local electric power source term to zero in the heat equation when *T*>100^{°}C. This allows us to approximate the material properties evolution linked to phase change and prevents very high non-physical temperatures.

### (e) Pressure equation and porous flow modelling

The porous blood flow used in the convection term of bio-heat equations is calculated from the gradient of the pressure, which is the solution of the Poisson equation. The blood flow, owing to the small vessel branches that are not visible as they are below the CT scan resolution, is modelled by assuming a porous flow boundary condition on the vessel wall surface. We can therefore approximate this porous flow with a pressure distribution on the vessel wall that, in turn, forms the boundary condition for the solution of the Poisson equation. Indeed, fixed positive (respectively, negative) pressure values on input (respectively, output) vessels walls allow the definition of a reasonable approximation of porous blood flow. This modelling choice has the advantage of being numerically stable as no point sources have to be introduced in the Poisson equation for the pressure. The values of pressure at boundaries are guided by the data of the regional hepatic blood flow *q*=26 ml 100 g^{−1} min^{−1}, which is then used to define an average target velocity across the liver tissue porous flow region. In addition to the fixed pressure on the vessel walls, a zero pressure flux boundary condition is applied on the liver outer wall.

### (f) Cell death model implementation: resolution of equations and modification of heat equation

The three-state tissue death model for necrosis is computed by solving the coupled transient equations expressed in alive and dead variables with initial states *A*=0.99 and *D*=0 over the computational domain. The vulnerable variable is then easily derived. A fixed-point method is implemented for nonlinearities arising from the temperature dependence of the forward rate constant in the mathematical model. The solution of this model is then used to modify some of the material parameters in the temperature equations. Indeed, to model the poor convection in carbonized cells, the convection term is locally set to zero when *D*>0.8, which is the critical threshold of cell viability in the model. A modified heat capacity value is also used locally below this threshold.

### (g) Parallel computing

A first parallel implementation of the finite-element code, based on mesh partition, uses the message passing interface standard for inter-processes communication, and allows the distribution of the resolution of different parts of the computational domain to multiple processors. The Elmer-integrated meshing tool (Elmer Grid), using the METIS library, creates the suitable mesh partitions. The use of the METIS library is well adapted to complicated geometries, as it ensures that the number of shared nodes between the partitions, and consequently also the amount of inter-processor communication during the parallel computation is minimized. An optimized parallel implementation is currently under development, one that takes advantage of potential model simplifications, guided by the more extensive validation of numerical results that are being gathered currently.

### (h) Meshing

The volume-meshing strategy for the complicated liver vessel tree structure is based on a linear four-noded tetrahedral element. The open source pre-processor EnGrid has been substantially modified to provide different surface mesh preparation tools for the segmented liver/vessel wall, which has a jagged and castellated surface mesh. The resulting surface mesh of the integrated portal vein, hepatic vein, hepatic artery and liver wall provides a watertight boundary for the volume meshing. The volume meshing is implemented via EnGrid with a manual mesh refinement. The adaptive mesh refinement is currently under development. The surface modification includes smoothing and coarsening of the segmented mesh in order to optimize the element count and the solution time. These surface modification tools are part of the simulation development process and although they will be available to the clinician, it is expected that we will have decided upon a suitable surface mesh density by the end of the image-based multi-scale physiological planning for ablation cancer treatment (IMPPACT) project. The surface mesh density will then be fixed and the final software will allow us to more easily automate the surface preparation and volume-meshing stages into the intervention planning system infrastructure. This is the only practical approach for clinicians who have no experience in the use of computational engineering simulation tools.

Typical surface mesh sizes are of the order of several hundred thousand triangular elements with resulting volume meshes of the order of 1–3 million tetrahedral elements. The mesh size is partly dependent on the number of small vessel branches near the ablation zone. The more of the small vessel branches there are, the more intricate the surface mesh and thereby the larger the number of elements in the volume mesh. The resolution from the CT scan, and thus the resolution of the segmentation, determines the vessel tree resolution. We expect the patient data to have fewer small branches than the pig data owing to the different contrast enhancement procedures used in clinical practice; therefore, the patient meshes and solution times ought to be faster and easier to implement than the pig simulation models.

We present here the results of the numerical simulation of two RFA processes on pig livers, namely RFA#26 and RFA#42. In both cases, the experimental data for electrical power source and geometry are available as inputs of the model. Indeed, the total electrical power over time is tracked by the RFA generator, and the segmentation of CT scans provides a detailed liver geometry (outer wall, hepatic artery, hepatic vein and portal vein). Moreover, the positions of the needle tips are also extracted. From these data, the heating source term and porous convection terms in the bio-heat equations can be defined, as well as the boundary conditions for the pressure solver.

Tetrahedral meshes were created here from the watertight surface mesh of the liver outer wall and inner vessel trees. After extraction of a refined mesh in a sphere around the segmented lesion, the resulting meshes for RFA#26 and RFA#42 contain 800 000 and 1 800 000 tetrahedral elements, respectively. These meshes describe the small vessel branches around the heating source in order to take into account their cooling effect in the RFA process. The density of the mesh is thus conditioned both by the geometry (small vessels branches inducing small surface elements) and by the numerical resolution of the model (electrical heating at smoothed point sources combined with porous flow convection). In particular, a fine mesh needs to be defined around the heat sources to eliminate numerical oscillations when porous flow convection is applied.

## 5. Results and discussion

In §§2–4, we have presented the necessary components of the tool: the mathematical model and its implementation, and the generation of the individual subject geometry. We now present the results of the simulations and compare them with the experimental measurements, enabling us to perform a quantitative validation of the model predictions.

The implementation of the mathematical model has been validated over several test cases provided by RFA experiments on pig livers. Control of the ablated tissue volume is the key feature of the patient-specific intervention planning system for use by the clinicians. The segmentation of the CT scan of the ablated tissue is therefore used for the comparison between the experimental lesion and the numerically predicted one. This segmentation gives the ablated tissue volume just after the RFA process, and one week/month later, in order to take into account the possible healing of cells. According to the mathematical model for cell death, the critical threshold of cell viability *D*>0.8 is then used to compute the boundary of the numerically predicted lesion.

The temperature distribution in the tissue, which is a function of the power input magnitude and distribution, the cooling effect of the vessels, the tissue and blood material properties and the tissue perfusion, can also be compared with the experimental data at some precise locations. In fact, five thermocouples placed at the ends of the RFA probe tips track the tissue temperature throughout the RFA process. This tracking of the temperature provides a control of the intermediate results/variables of the mathematical model, and is essential for the validation of the cell death model. More extensive temperature measurements would, however, be very valuable to provide a direct validation of the bio-heat model (as opposed to the whole model validation performed here) and this will be the subject of future work.

The computation is run in parallel on a cluster with three nodes, each node having two CPUs. Under these conditions, the computation of each RFA process (approx. 10 min in real time) needs approximately 4 h to complete. The speed-up factor with sequential computation is optimal for linear system assembly steps. However, the scalability of the linear solvers has still to be improved to increase the global speed-up factor. The results (temperatures, blood velocities, dead elements, etc.) are written out at specified time steps for visualization. The tissue temperature for RFA#26 at 600 s after ablation start is shown in figure 4*a* to illustrate both the localized nature of the heating and the large amount of cooling provided by the vasculature in particular regions of the liver. It is clear that the interaction between the heating and the temperature fields is highly complex, justifying the need for such a sophisticated model.

The remainder of figure 4 shows the results of both the predicted dead zone and the experimentally measured lesion size for figure 4*b* (RFA#26) and figure 4*c* (RFA#42). The sections of the vessel trees in the surrounding region are also shown, illustrating the importance of the heat sink effect. This cooling influence of the vessels modifies locally the temperature distribution and in turn reduces the predicted ablation zone. The correct thermal boundary conditions on the vessel walls near the ablation zone are therefore important for the accuracy of the final prediction. In particular, the application of convective boundary conditions for temperature allows the predicted ablation volume to wrap around some thermally non-significant vessels, as shown at the top of the figure for RFA#42.

These results show very good agreement for the size and shape of the ablated zone between the predicted results and the actual data. A more quantitative comparison is performed by computing the maximal deviation between the predicted and the actual lesion surfaces: likely to be the parameter of greatest clinical importance in ensuring that the predictions are sufficiently accurate for clinical practice. The values calculated here are 4.5 and 4.2 mm for the RFA#26 and RFA#42 pigs, respectively, with lesion characteristic lengths of 27.5 and 27 mm. Furthermore, the numerical temperatures at the tips ends are also in good agreement with the experimental data (data not shown). The volumes of the segmented and predicted lesions are given in table 2, again showing good agreement. Owing to the low thermal conductivity of the blood and the tissue, the temperature variation is also quite localized close to the heating source, which allows us to limit the computational domain to a volume approximately 35 mm radius from the centre of the ablation zone. This has a significant impact on the computational requirements of the simulations.

It should be noted that these results, although very encouraging, are still only preliminary, with only two comparisons being made. Now that the loop is completed, though, we will be performing further validation using further animal models. This will enable us to gain a greater insight into the likely quantitative error of our model. This dataset will also be invaluable in allowing us to refine the mathematical model to improve the accuracy of prediction. It is likely that there will be a lower bound of error, below which these models cannot go, however, owing to inherent uncertainties and errors in model parameter values and image quality. Determining this lower bound will be a very valuable contribution to this area, providing as it will be the inherent ability of mathematical modelling to contribute to clinical practice and ensuring that the clinicians are aware of the likely errors in predictions.

One of the greatest likely sources of error is the values of model parameters, such as thermal conductivity, which will vary from patient to patient and even spatially within the liver. The literature is very sparse on this, hence for a brief introduction to the limitations of the existing experimental data, e.g. the review by Payne *et al*. [3]. Considerable further work will need to be performed to reduce the errors associated with such variability, in particular, examining the variation in material properties with disease, such as cirrhosis, which is known to affect the flow field and is likely to strongly affect the material parameters. In addition, tumour flow is highly abnormal, often with stagnant flow; incorrect assumptions on the blood flow velocity field may thus overestimate the cooling effect significantly. Likewise, tumour hypoxia will be an additional factor in making the tumour response very different from the response of healthy liver tissue. As the tumour cells die, toxins are released, and these also have an effect on the remainder of the tumour.

Once we start to obtain clinical data from human patients, it will, however, be possible to examine the effects of such errors on the final prediction of lesion size; this sensitivity analysis will also be a valuable tool towards obtaining further experimental data where it is most important and thus needed. This further illustrates the importance of quantitative mathematical modelling and validation in helping to understand which biological questions most need answering and driving biological research forward in combination.

## 6. Conclusions

Comparison between model predictions and experimental data acquired on healthy animal models shows a very good agreement for the distribution of the ablated zone between the predicted and actual data. By carrying out a direct comparison between the two, it is possible to perform a quantitative validation of the model: a vital part of the development of any system. This is the first such attempt in the liver to close the gap between mathematical modelling and experimental validation of ablation treatment of tumours *in vivo*. As the project proceeds and more experimental data are gathered, further refinement of the model will be performed, in particular, to include animal-/patient-specific information, such as cirrhosis level. Validation in humans will then be performed as the project moves towards practical implementation, with substantial potential for widespread clinical use of a patient-specific intervention planning system. This tool will also provide quantitative results (particularly from validation) that will be very valuable in the development of similar systems for other body organs.

## Acknowledgements

The research leading to these results has received funding from the European Community's Seventh Framework Programme under grant agreement no. 223877, project IMPPACT.

## Footnotes

One contribution of 11 to a Theme Issue ‘Towards the virtual physiological human: mathematical and computational case studies’.

- This journal is © 2011 The Royal Society