Robust GPU-based Virtual Reality Simulation of Radio Frequency Ablations for Various Needle Geometries and Locations

TL;DR

This paper introduces a GPU-accelerated, real-time simulation method for radio-frequency ablations in liver cancer treatment, improving accuracy and safety in virtual reality surgical planning and guidance.

Contribution

It presents a novel GPU-based simulation approach that outperforms existing methods in accuracy, robustness, and real-time performance for radio-frequency ablation planning.

Findings

Achieves >480 Hz rendering on Nvidia 1080 Ti GPU

Shows significant correlation improvements over literature methods

Ensures no failure modes after initial 10 seconds

Abstract

Purpose: Radio-frequency ablations play an important role in the therapy of malignant liver lesions. The navigation of a needle to the lesion poses a challenge for both the trainees and intervening physicians. Methods: This publication presents a new GPU-based, accurate method for the simulation of radio-frequency ablations for lesions at the needle tip in general and for an existing visuo-haptic 4D VR simulator. The method is implemented real-time capable with Nvidia CUDA. Results: It performs better than a literature method concerning the theoretical characteristic of monotonic convergence of the bioheat PDE and a in vitro gold standard with significant improvements (p < 0.05) in terms of Pearson correlations. It shows no failure modes or theoretically inconsistent individual simulation results after the initial phase of 10 seconds. On the Nvidia 1080 Ti GPU it achieves a very high frame rendering performance of >480 Hz. Conclusion: Our method provides a more robust and safer real-time ablation planning and intraoperative guidance technique, especially avoiding the over-estimation of the ablated tissue death zone, which is risky for the patient in terms of tumor recurrence. Future in vitro measurements and optimization shall further improve the conservative estimate.

Tables

Table 1. Table 1 : Kernel-layout ( x ⋅ y ⋅ z ⋅ 𝑥 𝑦 𝑧 x\cdot y\cdot z -Threads/block) vs. runtime for 6000 images.

Kernel conf.	Time for 6000 images (s)	Time/image (s)	FPS (Hz)
$64\cdot 3\cdot 3$	10.8	0.0018	556
$256\cdot 3\cdot 1$	12.3	0.00205	488
$1024\cdot 1\cdot 1$	12.23	0.00205	490

Equations

$$\dot{T}=k{\nabla }^{2}T$$

$$T_f(x,t)=\frac{1}{\sqrt{4\pi kt}}\exp \left(-\frac{x^2}{4kt}\right).$$

$$\{\begin{array}{ll}\dot{T}=k\cdot \frac{\partial T}{\partial x^2}+f(x,t)& \text{additional terms f (AT)}\\ T(x,0)=g(x)& \text{initial condition (IC)}\\ T(0,t)=h(t)& \text{boundary condition (BC)}\end{array}$$

$$T_{BC}(x,t)={\int }_0^t\frac{x}{\sqrt{4\pi k(t-s{)}^3}}\exp \left(-\frac{x^2}{4k(t-s)}\right)h(s)ds,\forall x>0$$

$$T_0\cdot \int \frac{x}{\sqrt{4\pi k{t}^3}}\exp \left(-\frac{x^2}{4kt}\right)dt=-T_0\cdot \frac{\sqrt{k}{t}^{3/2}\cdot \text{erf}\left(\frac{x}{2\sqrt{k}}\right)}{\sqrt{k{t}^3}}+\text{const.}$$

$$T_{BC}(x,t)=T_0-T_0\cdot \frac{\sqrt{k}{t}^{3/2}\cdot \text{erf}\left(\frac{x}{2\sqrt{k}}\right)}{\sqrt{k{t}^3}}$$

$$\alpha (t)=\frac{1}{ln(c(0)/c(t))}={\left[A{\int }_0^t\text{exp}\left(-\frac{\Delta E}{RT(t)}\right)dt\right]}^{-1}$$

$$\rho c_p\frac{\partial T}{\partial t}=▽(K▽T)+\sigma |▽V{|}^2-{\rho }_b{w}_b{c}_b\cdot \alpha (t)(T-T_a).$$

$$\rho c_p\frac{\partial T}{\partial t}=▽(K▽T)+w_b{c}_b(T_a-T)+Q_m.$$

$$\rho c_p\frac{\partial T}{\partial t}=K\cdot \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)+w_b{c}_b(T_a-T)+Q_m.$$

$$\frac{{\partial }^2{T}_i}{\partial x^2}\approx \frac{T(x_{i+1})-2T(x_i)+T(x_{i-1})}{(△x{)}^2}.$$

$$\rho c_p\frac{T_{i,j,k}^{n+1}-T_{i,j,k}^n}{△t}=K\cdot \frac{T_{i-1,j,k}^n-2T_{i,j,k}^n+T_{i+1,j,k}^n}{(△x{)}^2}+\frac{T_{i,j-1,k}^n-2T_{i,j,k}^n+T_{i,j+1,k}^n}{(△y{)}^2}$$

$$+\frac{T_{i,j,k-1}^n-2T_{i,j,k}^n+T_{i,j,k+1}^n}{(△z{)}^2})+(w_{b_{i,j,k}}{c}_{b_{i,j,k}}(T_a^n-T_{i,j,k}^n)+Q_{m_{i,j,k}})$$

$$\frac{{\partial }^{2}T}{\partial x^2}=0,\text{for neighbouring voxel outside of liver (towards vessel or tip)}.$$

$$T\left({\mathbf{x}}_t\right)=T_{ref}\left({\mathbf{x}}_t+{\hat{u}}^{-1}({\mathbf{x}}_t,t)\right).$$

Full text

Methods: This publication presents a new GPU-based, accurate method for the simulation of radio-frequency ablations for lesions at the needle tip in general and for an existing visuo-haptic 4D VR simulator. The method is implemented real-time capable with Nvidia CUDA.

Results: It performs better than a literature method concerning the theoretical characteristic of monotonic convergence of the bioheat PDE and a in vitro gold standard with significant improvements ($p<0.05$) in terms of Pearson correlations. It shows no failure modes or theoretically inconsistent individual simulation results after the initial phase of 10 seconds. On the Nvidia 1080 Ti GPU it achieves a very high frame rendering performance of $>$480 Hz.

Conclusion: Our method provides a more robust and safer real-time ablation planning and intraoperative guidance technique, especially avoiding the overestimation of the ablated tissue death zone, which is risky for the patient in terms of tumor recurrence. Future in vitro measurements and optimization shall further improve the conservative estimate.

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 lower-risk radio-frequency ablation (RFA) 2675-01 is possible (Fig. 1(a)). RFA is the most common goldberg2002comparison and geometrically flexible ablation technique (see Fig. 1(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 $-40\text{\,}\mathrm{\SIUnitSymbolCelsius}$ (cryoablation) or exceeding $60\text{\,}\mathrm{\SIUnitSymbolCelsius}$ 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. 1(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 :

[TABLE]

describes the temperature $T$ distribution in a solid medium over time. $k$ is a constant describing the temperature (substance or energy) spread rate. Now, we can introduce a border condition with an initial point source $T(x,0)=g(x)=\delta(x)$ of heat at a known position, i.e. 0 with heat $1\text{\,}\mathrm{\SIUnitSymbolCelsius}$. 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:

[TABLE]

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 :

[TABLE]

where the first line is the fundamental Laplacian with other additive terms $f(x,t)$, IC corresponds to a initial temperature distribution (later: image) at time point zero and BC suits to represent a constant heat source $h(t)=T_{0}$. In this chapter for brevity, we set $g(x)=0$ 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 $f(x,t)=0,g(x)=0,h(t)=T_{0}$, 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 $T(0,t)=h(t)=\delta(x)\cdot$90\text{,}\mathrm{\SIUnitSymbolCelsius}$$. The solution of the time convolution:

[TABLE]

can be found using the indefinite integral - as $h(s)=T_{0}=const.$, i.e.:

[TABLE]

and finally yields with integration limits applied and $const.=T_{0}$:

[TABLE]

This characteristic analytic solution $T=T_{BC}$ is plotted with $h(t)=T_{0}=$90\text{,}\mathrm{\SIUnitSymbolCelsius}$$ and $k=0.52$ in Fig. 2(a) over space $x\geq 0$ and time $t\geq 0$. The solution shows monotonic convergence of temperature (Figs. 2(b), 2(c)).

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:

[TABLE]

with symbols $c(t)$ as the healthy cell concentration at time $t$, $R$ as the universal gas constant ($8.31\text{\,}\frac{\mathrm{J}}{\mathrm{mol}\text{\,}\mathrm{K}}$), $A$ as a frequency factor for the kinetic expression ($\mathrm{Hz}$), $\Delta E$ denotes the activation energy for the irreversible cell damage reaction ($\frac{\mathrm{J}}{\mathrm{mol}}$) and $T$ denotes the tissue temperature ($\mathrm{\SIUnitSymbolCelsius}$) 2675-05 ; 2675-20 : Damaged tissue with no cooling by blood flow corresponds to $\alpha=0$ at some later time point, and fully functional tissue is represented by $\alpha=1$ at $t=0$ as shown in Fig. 3. $\alpha$ 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 $37\text{\,}\mathrm{\SIUnitSymbolCelsius}$. 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 :

[TABLE]

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:

[TABLE]

The variables denote the resistive heat transfer at Linte’s needle tip $\sigma|\triangledown V|^{2}$ (with $V$ as the electric potential ($\mathrm{V}$)), the blood density $\rho_{p}$ ($\frac{\mathrm{kg}}{{\mathrm{m}}^{3}}$) and the tissue state coefficient $\alpha(t)$ (Sec. 2.1.2) from Eq. 8. Eq. 9 incorporates the metabolic heat generation $Q_{m}$ of liver cells ($10\,714\text{\,}\frac{\mathrm{W}}{{\mathrm{m}}^{3}}$) 2675-03 . In both equations, $\rho$ the specific density of liver tissue ($1079\text{\,}\frac{\mathrm{kg}}{{\mathrm{m}}^{3}}$), $c_{p}$ denotes the liver tissue heat capacity ($3540\text{\,}\frac{\mathrm{J}}{\mathrm{kg}\text{\,}\mathrm{\SIUnitSymbolCelsius}}$), $t$ the time ($\mathrm{s}$), $K$ the liver thermal conductivity ($0.52\text{\,}\frac{\mathrm{W}}{\mathrm{m}\text{\,}\mathrm{\SIUnitSymbolCelsius}}$), $T$ the temperature ($\mathrm{\SIUnitSymbolCelsius}$), $w_{b}$ the liver blood flow rate ($16.687\text{\,}\frac{\mathrm{kg}}{{\mathrm{m}}^{3}\text{\,}\mathrm{s}}$), $c_{b}$ the heat capacity of the blood ($3617\text{\,}\frac{\mathrm{J}}{\mathrm{kg}\text{\,}\mathrm{\SIUnitSymbolCelsius}}$), and $T_{a}$ the blood temperature ($37\text{\,}\mathrm{\SIUnitSymbolCelsius}$).

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 $Q_{m}$ and neglect external cooling $Q_{p}$ (air) ($\frac{\mathrm{W}}{{\mathrm{m}}^{3}}$) 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 $q_{p}$ in ($\frac{\mathrm{W}}{{\mathrm{m}}^{3}}$), 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 $Q_{m}$ 2675-03 . Finally, Linte’s tissue decay term from Eq. 7 contains a singularity regarding $c(0)$ and yields high values near $t=0$.

For 3D continuous space, the partial differential equation can be written as:

[TABLE]

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 $x$:

[TABLE]

$\vartriangle x$ denotes the voxel distance between two voxels in the x-direction in millimeter. Neighbours to a voxel with index $i$ are the predecessor $i-1$ and the successor $i+1$ in the x-direction. The voxel distances $\vartriangle y$ and $\vartriangle z$ as well as indices $j$ OliveGreenand $k$ are used accordingly in the equations. Now, Eq. 10 can be transferred by substitution into 3D voxel image space $T^{n}_{i,j,k}$ for implementation:

[TABLE]

Inside a temperature image $T$, $T^{n}_{i,j,k}$ addresses the temperature value stored in a voxels at the position $i,j,k$ in image coordinates at the time iteration index $n\in N_{0}$. $w_{b_{i,j,k}}$, $c_{b_{i,j,k}}$ and $Q_{m_{i,j,k}}$ 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 ($37\text{\,}\mathrm{\SIUnitSymbolCelsius}$). 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:

[TABLE]

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 $\mathrm{X}$ with a diffeomorphic motion vector field $\hat{u}$ into the breathing motion’s curved coordinate space $\mathrm{x}$ (Fig. 4):

[TABLE]

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 $u_{phase}$. In a keyframe approach, consecutive vector fields $u_{phase}$ and $u_{phase+1}$ 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 $\hat{u}(\mathbf{x}_{t},t)$, 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 2675-21 .

The used subsampled and segmented CT data set (dimensions 2562x234, spacing $1.552\text{\,}\mathrm{mm}$2 x$2\text{\,}\mathrm{mm}$) is converted into a first thermal image ($t=n=0$) 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 $0.5\text{\,}{\mathrm{mm}}^{3}$ voxels. The $2.5\text{\,}\mathrm{mm}$ 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 $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and $60\text{\,}\mathrm{\SIUnitSymbolCelsius}$ are applied. Simulation time steps of $0.1\text{\,}\mathrm{s}$ 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 ($60\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$) 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 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 $5\text{\,}\mathrm{mm}$ and an ablation heat of $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$ 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 in vitro (gold standard) and our simulation did not exceed the tissue death threshold of $42.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ (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 $42.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$. The model of Linte 2675-05 is on average farther away from the gold standard with a larger range of fluctuation of SD= $0.99\text{\,}\mathrm{\SIUnitSymbolCelsius}$. 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 $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$ or $60\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and $5\text{\,}\mathrm{mm}$ radius (Fig. 8), Linte critically or totally fails. Overall, our model results are similar to the in vitro temperatures, which reflects in significantly high ($p<0.005$) Pearson correlation coefficients (Fig. 13).

The simulation of the model of this work proves a theoretically sound monotonically converging heating process (Figs. 14(a), 14(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 (r>$$5\text{\,}\mathrm{mm}). 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 $42.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ and red denotes $90\text{\,}\mathrm{\SIUnitSymbolCelsius}$.

3.2 Rendering-Performance

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 $12\%$ faster runtime of the configuration with $64\cdot 3\cdot 3$ 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 2675-05 . With smaller curvature errors, the model presented here as well as the in vitro determined temperatures are below the threshold of $42.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ (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 (FPS$>$24 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 $Q_{m}$ 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 $\alpha$-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 ($t=0$), what exactly is the formula for the function $c(t)$, 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 $c(t)$ 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).

5 Figures

6 Tables

7 Algorithms