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 lower-risk radio-frequency ablation (RFA) 2675-01 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 2675-10 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 2675-24 . Radiation tip brachytherapy is more common in the lower abdomen targeting the cervix or prostate 2675-15 ; 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 2675-25 . 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 radio-freqency generator heat up the tissue around the needle tip. Within the 4D VR needle ablation simulator from our group 2675-08 ; 2675-09 ; 2675-10 ; 2675-14 ; 2675-16 ; 2675-17 , we use an efficiently parallel calculable temperature model 2675-02 specialized to liver lesions. The new GPU-based algorithm also compares qualitatively (patient safety) favorable to the literature 2675-05 ; 2675-18 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 (under-treatment) and the 100% reproducability of simulations (no statistical errors). Currently, the training of ablations is carried out in a risky set-up directly on the living patient. However, VR training and planning simulations such as ours 2675-13 ; 2675-08 are on the rise and safe over-real-time 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 :
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.1.1 General Formal Framework
Generally, a linear combination of solutions of other well-known 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 :
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:
can be found using the indefinite integral - as , i.e.:
and finally yields with integration limits applied and :
2.1.2 Cell Decay Terms from the Comparison Work
Linte et al. 2675-05 proposed a tissue state decay term based on an inverted Arrhenius equation:
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 () 2675-05 ; 2675-20 : 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 under-weighting the blood cooling.
2.2 Specialized Methods for Liver Ablation Simulation
The segmented regions of the virtual bodies 2675-11 ; 2675-06 ; 2675-07 are assigned a structure-specific 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 2675-02 . (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. 2675-05 ; 2675-02 :
Assessment of the relevant contribution of different terms of Eq. 8 vs. in vitro experiments 2675-05 has lead us to a more robust simplified simulation equation used in this work for liver tumor ablation:
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 () 2675-03 . 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 2675-02 . The simulation of breathing motion is optionally possible, see Sec. 2.3. Compared to Pennes 2675-02 , 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 liver-specific 2675-03 . 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:
For 3D continuous space, the partial differential equation can be written as:
The finite-difference method (FDM) with symmetric partial difference quotients is used to discretize the partial derivatives on a regular grid of voxels, i.e. here for :
denotes the voxel distance between two voxels in the x-direction in millimeter. Neighbours to a voxel with index are the predecessor and the successor in the x-direction. 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:
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 voxel-wise.
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:
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 2675-08 ; 2675-09 . The simulated temperature images are deformed in the reference coordinate space with a diffeomorphic motion vector field
with a diffeomorphic motion vector fieldinto the breathing motion’s curved coordinate space (Fig. 4):
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. Non-linear registrations are calculated between the phases resulting in inter-phase 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
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 .
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 2675-21 .
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 2675-03 .
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. 2675-05 ; 2675-18 . 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. 2675-05 as a voxel hemisphere by five voxels in xy-direction and three voxels in z-direction. Temperatures of and are applied. Simulation time steps of were selected.
Linte et al. 2675-05 ; 2675-18 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 gold-standard 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 ) 2675-18 . 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.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 under-planning 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
is a key experiment, where Linte’s method shows significant under-planning 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 measuredin 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 2675-05 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 2675-05 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).
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 minute111https://goo.gl/6VPxSU. Second in the context of the 4D-VR 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 222https://bit.ly/2HamKtO. In the movie footage, the yellow border of the ablation zone corresponds to and red denotes .
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 real-time 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.
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 2675-05 . 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 over-fulfilling frame rate (FPS24 Hz, x23) for our VR planning simulator 2675-13 and in general. In contrast to Linte et al. 2675-05 , 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 2675-05 won, but lost in the other three challenges for temperature (ablation zone) overestimation (Figs. 8, 9) or flat-line (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 re-implement 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 in-vitro experiments with significant sample size are planned by us without necessarily comparing to Linte et al.
Acknowledgements.Funding: DFG: MA 6791/1-1; Nvidia GPU Grant 2018 (Mastmeyer).
- (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 meta-analysis. 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.: Gpu-based visualization of deformable volumetric soft-tissue for real-time 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 image-based multiproxy palpation algorithm for patient-specific vr-simulation. 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 visuo-haptic 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 visuo-haptic 4D volume rendering using respiratory motion models. IEEE Transactions on Haptics 8(4), 371–383 (2015)
- (9) Goldberg, S.N.: Comparison of techniques for image-guided 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 radio-frequency ablation lesions for image-guided 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 image-guided arrhythmia therapy: A preliminary ex vivo demonstration. In: Augmented Environments for Computer-Assisted 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 visuo-haptic VR simulation using a generic patient atlas. Computer Methods and Programs in Biomedicine 132, 161–175 (2016)
Mastmeyer, A., Fortmeier, D., Handels, H.: Random forest classification of large volume structures for visuo-haptic 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.: Patch-based label fusion using local confidence-measures 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.: Model-based catheter segmentation in mri-images. In: International Conference on Medical Image Computing and Computer-Assisted Intervention – MICCAI (2015)
- (22) Mastmeyer, A., Pernelle, G., Ma, R., Barber, L., Kapur, T.: Accurate model-based segmentation of gynecologic brachytherapy catheter collections in mri-images. Medical image analysis 42, 173–188 (2017)
- (23) Mastmeyer, A., Wilms, M., Fortmeier, D., Schröder, J., Handels, H.: Real-time ultrasound simulation for training of US-guided 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.: Population-based 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(11-12), 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)
|Kernel conf.||Time for 6000 images (s)||Time/image (s)||FPS (Hz)|