1 Introduction
In liver cancer therapy besides the conservatively used open laparotomy goldberg2002comparison ; oshowo2003comparison , there are several options for the treatment of primary liver tumors.
The oral endoscopic access is chosen for lesion in or at the main ducts (e.g. CHD). The minimally invasive laparoscopic access (endoscopic keyhole surgery) is favorable for peripheral liver boundary tumors. For smaller tumors in the vicinity of hepatic arteries, the minimally invasive and lowerrisk radiofrequency ablation (RFA) 267501 is possible (Fig. 0(a)). RFA is the most common goldberg2002comparison and geometrically flexible ablation technique (see Fig. 0(b)) for focal tumors and therefore is the basis of this work. A needle containing an electrode is used to access the target, and the electrode is used to emit high frequency currents to the needle tip. Fluoroscopy or ultrasound 267510 image guidance helps to navigate the tip to the target lesion. Compared to open laparotomy, the survival chance is only 2% worse oshowo2003comparison . The long chained proteins of the cells disintegrate leading to the tissue ablation. Clinically, total necrosis instantaneously occurs at temperatures below (cryoablation) or exceeding for most cell types nikfarjam2005mechanisms .
Uncommonly used, microwave ablation is a rather new technique which heats tissue around a coil in the tip to remove cancer tissue 267524 . Radiation tip brachytherapy is more common in the lower abdomen targeting the cervix or prostate 267515 ; mastmeyer2015model . The treatment of spherical or ellipsoidal tumours is faster, subject to shorter fluoroscopy exposition and more painless chen2015efficacy with cryotherapy kuck2016cryoballoon , which destroys tumour cells with ice crystals from quick freezing and thawing 267525 . Transarterial chemoembolization (TACE) is another treatment option (yet more side effects) for rather small or otherwise unresectable tumours extending the lifetime of patients llovet2003systematic . Laser ablation based on dispersed laser light from fibers is another safe and successful thermal ablation technique izzo2003other .
The ablation principle is shown in Fig. 0(a). Currents from a radiofreqency generator heat up the tissue around the needle tip. Within the 4D VR needle ablation simulator from our group 267508 ; 267509 ; 267510 ; 267514 ; 267516 ; 267517 , we use an efficiently parallel calculable temperature model 267502 specialized to liver lesions. The new GPUbased algorithm also compares qualitatively (patient safety) favorable to the literature 267505 ; 267518 and theoretical characteristics. Keep in mind throughout the rest of the paper in RFA planning, overestimation of the heat zone can be deadly by finally leaving residual tumor cells (undertreatment) and the 100% reproducability of simulations (no statistical errors). Currently, the training of ablations is carried out in a risky setup directly on the living patient. However, VR training and planning simulations such as ours 267513 ; 267508 are on the rise and safe overrealtime ablation planning is highly relevant for intraoperative planning systems.
2 Material and Methods
The calculation of the heat distribution is based on virtual body models, which essentially consist of segmented 3D CT image data. To better found our argumentation, we first show some theoretical thoughts.
2.1 Theoretical Background
The fundamental solution of the second order hyperbolic partial differential heat equation (PDE) evans2012numerical :
(1) 
describes the temperature distribution in a solid medium over time. is a constant describing the temperature (substance or energy) spread rate. Now, we can introduce a border condition with an initial point source of heat at a known position, i.e. 0 with heat . This is known as the fundamental heat distribution problem solved by convolution.
The fundamental solution to the above problem is a typical Green’s function:
(2) 
2.1.1 General Formal Framework
Generally, a linear combination of solutions of other wellknown fundamental problem solutions subject to other elementary initial (IC, initial temperature distribution) or border conditions (BC, constant heat source) or slight modifications in the PDE term in Eq. 2 on the right hand side (AT) evans2012numerical :
(3) 
where the first line is the fundamental Laplacian with other additive terms , IC corresponds to a initial temperature distribution (later: image) at time point zero and BC suits to represent a constant heat source . In this chapter for brevity, we set and concentrate to solve the AT and BC problems for our needs.
For elaborating just the characteristics, we sum up the math: The solution with respect to the special terms, set here to , is a linear combination of appropriate Green’s function solutions resulting from the linear operation of convolutions of the type shown in Eq. 2 (Bell curves).
This translates perfect to the characteristic simulation problem of needle ablation, e.g. using a needle tip heat source . The solution of the time convolution:
(4) 
can be found using the indefinite integral  as , i.e.:
(5) 
and finally yields with integration limits applied and :
(6) 
2.1.2 Cell Decay Terms from the Comparison Work
Linte et al. 267505 proposed a tissue state decay term based on an inverted Arrhenius equation:
(7) 
with symbols as the healthy cell concentration at time , as the universal gas constant (), as a frequency factor for the kinetic expression (), denotes the activation energy for the irreversible cell damage reaction () and denotes the tissue temperature () 267505 ; 267520 : Damaged tissue with no cooling by blood flow corresponds to at some later time point, and fully functional tissue is represented by at as shown in Fig. 3. is a monotonically decreasing term from 1 to 0 underweighting the blood cooling.
2.2 Specialized Methods for Liver Ablation Simulation
The segmented regions of the virtual bodies 267511 ; 267506 ; 267507 are assigned a structurespecific temperature, e.g. a liver blood vessel has a constant temperature of . The temperature flow in the human body consists of three constituents benzinger1969heat : (1) Heat conduction from warm to cold tissue. (2) Cooling by blood flow 267502 . (3) Metabolic heat generation. The model to compare against for temperature propagation using liver ablations in this work is the modified Bioheat equation as stated in Linte et al. 267505 ; 267502 :
(8) 
Assessment of the relevant contribution of different terms of Eq. 8 vs. in vitro experiments 267505 has lead us to a more robust simplified simulation equation used in this work for liver tumor ablation:
(9) 
The variables denote the resistive heat transfer at Linte’s needle tip (with as the electric potential ()), the blood density () and the tissue state coefficient (Sec. 2.1.2) from Eq. 8. Eq. 9 incorporates the metabolic heat generation of liver cells () 267503 . In both equations, the specific density of liver tissue (), denotes the liver tissue heat capacity (), the time (), the liver thermal conductivity (), the temperature (), the liver blood flow rate (), the heat capacity of the blood (), and the blood temperature ().
In contrast to Linte, in our planning system we model the needle tip or heat source with complex geometries (e.g. umbrellas), blood vessels as boundary conditions and breathing motion simulation, use a metabolic heat source and neglect external cooling (air) () from the historical Pennes formula 267502 . The simulation of breathing motion is optionally possible, see Sec. 2.3. Compared to Pennes 267502 , this work ignores the external sources spatial heat up phase in (), because in our situation in the liver there is no constant external cooling (air) and the blood cooling is already included in the middle term on the right side of the PDE. Moreover vs. Linte, we use a liverspecific 267503 . Finally, Linte’s tissue decay term from Eq. 7 contains a singularity regarding and yields high values near .
For 3D continuous space, the partial differential equation can be written as:
(10) 
The finitedifference method (FDM) with symmetric partial difference quotients is used to discretize the partial derivatives on a regular grid of voxels, i.e. here for :
(11) 
denotes the voxel distance between two voxels in the xdirection in millimeter. Neighbours to a voxel with index are the predecessor and the successor in the xdirection. The voxel distances and as well as indices OliveGreenand are used accordingly in the equations. Now, Eq. 10 can be transferred by substitution into 3D voxel image space for implementation:
(12)  
Inside a temperature image , addresses the temperature value stored in a voxels at the position in image coordinates at the time iteration index . , and denote the parameters declared for Eq. 10 voxelwise.
The Neumann boundary conditions for the blood vessel and needle tip voxels are ideally assumed as completely isolating, i.e. maintaining a constant temperature (). This is based on the reasonable assumption that any heating of the blood is immediately carried away by blood flow or constant tip heat by external energy:
(13) 
2.3 Consideration of Breathing Motion
In the 4D simulation of a liver puncture, the movement of the thorax, and here mainly the diaphragm and upper abdomen, during breathing is taken into account 267508 ; 267509 . The simulated temperature images are deformed in the reference coordinate space
with a diffeomorphic motion vector field
into the breathing motion’s curved coordinate space (Fig. 4):(14) 
The motion vector fields are computed from consecutive 3D CT phase images (static breathing states) taken from a breathing cycle scanned over time in a 4D CT data set. Nonlinear registrations are calculated between the phases resulting in interphase compensation motion vector fields . In a keyframe approach, consecutive vector fields and
serve to interpolate the motion of the cycle as scanned. In a linear regression approach, the vector fields (regressand) and a breathing surrogate signal (regressor, e.g. breathing volume) yield a time variant vector field
, optionally with surrogate induced irregular breathing motion mastmeyer2017interpatient ; mastmeyer2018population .2.4 Implementation
The Bioheat equation is a computationally intensive formula. Thus, the Bioheat equation is implemented using Nvidia CUDA and runs massively parallelized on a GPU with thousands of computing cores 267521 .
The used subsampled and segmented CT data set (dimensions 256x234, spacing x) is converted into a first thermal image () for Alg. 1. The voxels’ label values of the segmented data set are replaced by the temperatures of the denoted tissue 267503 .
Upon this start image (Fig. 5), the heat equation Eq. 12 is calculated and yields a temperature update value per voxel. In our VR simulator, an image double buffer priem1996apparatus is used to avoid rendering artifacts from unfinished calculations. Double buffering switches back and forth between two buffers using memory pointers, see Alg. 2.
2.5 Evaluation Methods
The temperature propagation simulation proposed here has been compared with in vitro measurements on real tissue and the simulation by Linte et al. 267505 ; 267518 . To compare in a fair manner, the model of this work is calculated on the same high resolution voxels. The geometry of the needle tip is modelled congruent to Linte et al. 267505 as a voxel hemisphere by five voxels in xydirection and three voxels in zdirection. Temperatures of and are applied. Simulation time steps of were selected.
Linte et al. 267505 ; 267518 conducted clinically relevant in vitro experiments. Alas, they were not reproducible by us to create more data samples, mainly either because the algorithm was not reproducable or made available from the authors work. The goldstandard for the comparison experiments in Linte’s in vitro experiments results from repeated measurements using fresh tissue samples. Symmetrically arranged temperature probes were placed in 2.5 mm and 5 mm distance from the heat source. Four according reference mean curves result for two distances vs. two ablation temperatures ( and ) 267518 . Our according in silico simulation measurement points are shown in Fig. 6.
The ”thread per block configuration” performance of the method is measured in frames per second (FPS). For measurements, a consumer Nvidia GTX 1080 GPU is used with 8 Gb memory and 2560 CUDA cores.
3 Results
First the quantitative results (Figs. 7, 8, 9, 10), then qualitative simulation examples near a blood vessel and using a complex ablation electrode geometry are shown (Figs. 11, 12).
3.1 Simulation of Heat Propagation
The experiment with a radius of and an ablation heat of
is a key experiment, where Linte’s method shows significant underplanning in terms of ablation time (vice versa: ablation zone predicted too big). His method is unsafe and would harm the patient. It is important to note, simulation results are absolutely reproducible without standard deviations. The temperatures measured
in vitro (gold standard) and our simulation did not exceed the tissue death threshold of (Fig. 8) after one minute. With smaller errors the model presented here as well as the in vitro determined temperatures are below the threshold of . The model of Linte 267505 is on average farther away from the gold standard with a larger range of fluctuation of SD= . The model from this paper delineates the temperature distribution for dying cells more accurately, robustly and safer than the comparative model 267505 as no failures (Fig. 10) or overestimations are observed (Figs. 8, 9). Statistically, this work’s model shows a systematically more similar, stable and theretically sound temperature curve to the gold standard. Summing up, in the key experiments with or and radius (Fig. 8), Linte critically or totally fails. Overall, our model results are similar to the in vitro temperatures, which reflects in significantly high () Pearson correlation coefficients (Fig. 13).The simulation of the model of this work proves a theoretically sound monotonically converging heating process (Figs. 13(a), 13(b)).
In the qualitative results and without the influence of a tissue border, the propagation of the cell death zone converges to a spherical shape ((Fig. 11a)). If another colder tissue structure influences the heat zone, its shape adapts to the boundary of the other tissue due to enhanced cooling (Fig. 11b, c). Finally, a consideration of the simulation of complex needle geometries such as umbrella needle shapes (Fig. 12) is interesting for the treatment of larger tumors (). First the tissue between the individual hot wires heats up in vitro and then the heat zone converges to a spherical shape. The first supplied movie demonstrates the development of the ablation zone over one minute^{1}^{1}1https://goo.gl/6VPxSU. Second in the context of the 4DVR simulator with simulated breathing motion according to Eq. 14, a tumor lesion with a yellow fringe is shown in the top left simulated fluoroscopy view ^{2}^{2}2https://bit.ly/2HamKtO. In the movie footage, the yellow border of the ablation zone corresponds to and red denotes .
3.2 RenderingPerformance
All kernel configurations (Tab. 1) achieve over 24 frames per second (FPS) by far, making it possible to achieve a smooth display needed in moving images display. All configurations could be used for realtime simulation of temperature propagation. The approximately faster runtime of the configuration with threads per block is due to the better use of the shared memory of CUDA in the calculation of the finite differences.
4 Discussion
A highly efficient, accurate and safe planning method for ablation prediction is proposed in this work. It has been shown that the model of this work achieves qualitatively plausible and quantitatively robust and safe results with different ablation temperatures and in terms of heat distribution. We take into account the influence of tissue boundaries (e.g. blood vessel cooling) and complex needle tip geometries. Unfortunately, we do neither have access to lab equipment or more in vitro/in silico data for the ablations nor an implementation of Linte’s algorithm. Summing up the quantitative evaluation, the model from this work predicts the temperature distribution for dying cells more accurately, robustly and reliably than the comparative model 267505 . With smaller curvature errors, the model presented here as well as the in vitro determined temperatures are below the threshold of (Fig. 8) in a key experiment and show a systematically more similar temperature curve vs. the gold standard and most importantly vs. theoretical curves. In the studies, the model of this work shows a conservative underestimation of temperature, but in contrast to the comparison work it shows no failures, overestimates and irregular fluctuations around the gold standard while being very near to our elaborated theoretical characteristics.
The presented model simulated single frames with a far overfulfilling frame rate (FPS24 Hz, x23) for our VR planning simulator 267513 and in general. In contrast to Linte et al. 267505 , the tissue state coefficient and thermal resistance are not considered, but the metabolic heat generation of the liver is. The heat zones from this work are consistently closer to reality, i.e. the radius of the heat zone is more accurate and safe preventing tumor recurrence. The systematic underestimation in Fig. 7 means a more conservative planning and is safer than an overestimation, which is critical to the patient (tumor recurrence). In only this experiment, the comparison method 267505 won, but lost in the other three challenges for temperature (ablation zone) overestimation (Figs. 8, 9) or flatline (Fig. 10). Note again, simulation result differences are absolutely reproducible with no standard deviations, i.e. they are highly significant per se  even if the mean gold standard curve is with statistical deviations.
The theoretical monotonically convergence characteristic is hurt in Linte’s method. In RFA planning, overestimation of the heat zone can be deadly by leaving residual tumor cells. We attribute this advantage of our method to the dubious term used by Linte depriving blood cooling too early. Additionally, several questions are not answered in his work: How is the cell concentration determined at the beginning of the ablation (), what exactly is the formula for the function , and how is the singularity avoided. These open issues are the strongest argument to use our method, which is transparently presented here. We were not able to get the source code from Linte or reimplement his method due to the imprecise method descriptions. We were left to compare against his provided experimental and simulation results. In the light of absolutely reproducible mathematical simulations, our method is a significantly better option. The mean comparison gold standard from Linte’s experiments results from multiple experiments. Our future research is the clarification what factors are represented in the experimental (linear) error correction terms to make our model even better and later come up with a fully theoretical approach once the remaining error factors are identified  e.g. we could use a better tissue decay term with a clear estimator for healthy cell concentration. In future research, new invitro experiments with significant sample size are planned by us without necessarily comparing to Linte et al.
Acknowledgements.
Funding: DFG: MA 6791/11; Nvidia GPU Grant 2018 (Mastmeyer).References
 (1) Benzinger, T.H.: Heat regulation: homeostasis of central temperature in man. Physiological reviews 49(4), 671–759 (1969)
 (2) Chen, Y.H., Lin, H., Xie, C.L., Zhang, X.T., Li, Y.G.: Efficacy comparison between cryoablation and radiofrequency ablation for patients with cavotricuspid valve isthmus dependent atrial flutter: a metaanalysis. Scientific reports 5, 10910 (2015)
 (3) Curley, S.: Radiofrequency ablation of malignant liver tumors. Annals of Surgical Oncology 10(4), 338–347 (2003)
 (4) Evans, G., Blackledge, J., Yardley, P.: Numerical methods for partial differential equations. Springer Science & Business Media (2012)
 (5) Fortmeier, D., Mastmeyer, A., Handels, H.: Gpubased visualization of deformable volumetric softtissue for realtime simulation of haptic needle insertion. In: Bildverarbeitung für die Medizin 2012, pp. 117–122. Springer, Berlin, Heidelberg (2012)
 (6) Fortmeier, D., Mastmeyer, A., Handels, H.: An imagebased multiproxy palpation algorithm for patientspecific vrsimulation. Studies in health technology and informatics pp. 107–113 (2014)
 (7) Fortmeier, D., Mastmeyer, A., Schröder, J., Handels, H.: A virtual reality system for PTCD simulation using direct visuohaptic rendering of partially segmented image data. IEEE Journal of Biomedical and Health Informatics 20(1), 355–366 (2016)
 (8) Fortmeier, D., Wilms, M., Mastmeyer, A., Handels, H.: Direct visuohaptic 4D volume rendering using respiratory motion models. IEEE Transactions on Haptics 8(4), 371–383 (2015)
 (9) Goldberg, S.N.: Comparison of techniques for imageguided ablation of focal liver tumors. Radiology 223(2), 304–307 (2002)
 (10) Hasgall, P., Di Gennaro, F., Baumgartner, C., Neufeld, E., Lloyd, B., Gosselin, M., Payne, D., Klingenböck, A., Kuster, N.: IT’IS Database for thermal and electromagnetic parameters of biological tissues (2018)
 (11) Izzo, F.: Other thermal ablation techniques: microwave and interstitial laser ablation of liver tumors. Annals of surgical oncology 10(5), 491–497 (2003)
 (12) Kuck, K.H., Brugada, J., Fürnkranz, A., Metzner, A., Ouyang, F., Chun, K.J., Elvan, A., Arentz, T., Bestehorn, K., Pocock, S.J., Albenque, J.P., Tondo, C.: Cryoballoon or radiofrequency ablation for paroxysmal atrial fibrillation. New England Journal of Medicine 374(23), 2235–2245 (2016)
 (13) Linte, C., Camp, J., Holmes, D., Rettmann, M., Packer, D., RA, R.: Toward modeling of radiofrequency ablation lesions for imageguided left atrial fibrillation therapy: model formulation and preliminary evaluation. Studies in health technology and informatics 184(11), 261–267 (2013)
 (14) Linte, C.A., Camp, J.J., Holmes, D.R., Rettmann, M.E., Robb Richard A.”, e.C.A., Chen, E.C.S., Berger, M.O., Moore, J.T., Holmes, D.R.: Modeling of radiofrequency ablation lesions for imageguided arrhythmia therapy: A preliminary ex vivo demonstration. In: Augmented Environments for ComputerAssisted Interventions, pp. 22–33. Springer Berlin Heidelberg, Berlin, Heidelberg (2013)
 (15) Llovet, J.M., Bruix, J.: Systematic review of randomized trials for unresectable hepatocellular carcinoma: chemoembolization improves survival. Hepatology 37(2), 429–442 (2003)
 (16) Mastmeyer, A., Fortmeier, D., Handels, H.: Direct haptic volume rendering in lumbar puncture simulation. In: Studies in Health Technology and Informatics: Medicine Meets Virtual Reality 19  MMVR 2012, Studies in Health Technology and Informatics, vol. 173, pp. 280–286. IOS Press (2012)
 (17) Mastmeyer, A., Fortmeier, D., Handels, H.: Efficient patient modeling for visuohaptic VR simulation using a generic patient atlas. Computer Methods and Programs in Biomedicine 132, 161–175 (2016)

(18)
Mastmeyer, A., Fortmeier, D., Handels, H.: Random forest classification of large volume structures for visuohaptic rendering in CT images.
In: Proc. SPIE Medical Imaging: Image Processing, vol. 9784, pp. 97842H–1–8. International Society for Optics and Photonics (2016)  (19) Mastmeyer, A., Fortmeier, D., Handels, H.: Evaluation of direct haptic 4d volume rendering of partially segmented data for liver puncture simulation. Nature Scientific Reports 7(1), 1–15 (2017)
 (20) Mastmeyer, A., Fortmeier, D., Maghsoudi, E., Simon, M., Handels, H.: Patchbased label fusion using local confidencemeasures and weak segmentations. In: Proc. SPIE Medical Imaging: Image Processing, pp. 86691N–1–11. International Society for Optics and Photonics, Orlando, USA (2013)
 (21) Mastmeyer, A., Pernelle Guillaume, B.L., Pieper, S., Fortmeier, D., Wells, S., Handels, H., Kapur, T.: Modelbased catheter segmentation in mriimages. In: International Conference on Medical Image Computing and ComputerAssisted Intervention – MICCAI (2015)
 (22) Mastmeyer, A., Pernelle, G., Ma, R., Barber, L., Kapur, T.: Accurate modelbased segmentation of gynecologic brachytherapy catheter collections in mriimages. Medical image analysis 42, 173–188 (2017)
 (23) Mastmeyer, A., Wilms, M., Fortmeier, D., Schröder, J., Handels, H.: Realtime ultrasound simulation for training of USguided needle insertion in breathing virtual patients. In: Studies in Health Technology and Informatics: Medicine Meets Virtual Reality 22  MMVR 2016, Studies in Health Technology and Informatics, vol. 220, pp. 219–226. IOS Press (2016)
 (24) Mastmeyer, A., Wilms, M., Handels, H.: Interpatient respiratory motion model transfer for virtual reality simulations of liver punctures. Journal of World Society of Computer Graphics  WSCG 25(1), 1–10 (2017)
 (25) Mastmeyer, A., Wilms, M., Handels, H.: Populationbased respiratory 4d motion atlas construction and its application for vr simulations of liver punctures. In: SPIE Medical Imaging 2018: Image Processing, vol. 10574, p. 1057417. International Society for Optics and Photonics (2018)
 (26) Meloni, M.F., Chiang, J., Laeseke, P.F., Dietrich, C.F., Sannino, A., Solbiati, M., Nocerino, E., Brace, C.L., Lee Jr, F.T.: Microwave ablation in primary and secondary liver tumours: technical and clinical approaches. International Journal of Hyperthermia 33(1), 15–24 (2017)
 (27) Nickolls, J., Buck, I., Garland, M., Skadron, K.: Scalable parallel programming with cuda. Queue 6(2), 40–53 (2008)
 (28) Nikfarjam, M., Muralidharan, V., Christophi, C.: Mechanisms of focal heat destruction of liver tumors. Journal of Surgical Research 127(2), 208–223 (2005)
 (29) Niu, L.Z., Li, J.L., Xu, K.C.: Percutaneous cryoablation for liver cancer. Journal of clinical and translational hepatology 2(3), 182 (2014)
 (30) Oshowo, A., Gillams, A., Harrison, E., Lees, W., Taylor, I.: Comparison of resection and radiofrequency ablation for treatment of solitary colorectal liver metastases. British journal of surgery 90(10), 1240–1243 (2003)
 (31) Pennes, H.: Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of Applied Physiology 1(2), 93–122 (1948)
 (32) Priem, C., Malachowsky, C., McIntyre, B., Moffat, G.: Apparatus for selecting frame buffers for display in a double buffered display system (1996). US Patent 5,543,824
 (33) Shen, W., Zhang, J., Yang, F.: Modeling and numerical simulation of bioheat transfer and biomechanics in soft tissue. Mathematical and Computer Modelling 41(1112), 1251–1265 (2005)
 (34) Werner, J., Buse, M.: Temperature profiles with respect to inhomogeneity and geometry of the human body. Journal of Applied Physiology 65(3), 1110–1118 (1988)
5 Figures
6 Tables
Kernel conf.  Time for 6000 images (s)  Time/image (s)  FPS (Hz) 

10.8  0.0018  556  
12.3  0.00205  488  
12.23  0.00205  490 
Comments
There are no comments yet.