1 Introduction

Cryoablation is an alternative to open surgery in the treatment of abdominal tumors. It consists in destroying the malignant cells using an extreme cold propagated from the tip of multiple cryoprobes inserted in the abdomen of the patient. This kind of interventions being lighter than open surgery and providing a faster recovery, it has become popular in the past decades to treat tumors of a reasonable size. However, a strong expertise is necessary to perform it, due to the difficulty to find a safe and optimal position for the cryoprobes, and to predict accurately the final shape of the iceball produced by the combination of the multiple freezing sources. This paper focuses on the second issue, with the aim of assisting the surgeon in predicting the final iceball using computer simulation.

Early approaches for computer assistance to cryoablation planning were mostly using ellipsoids to represent the shape of the iceball, with the idea of mimicking the theoretical measurements provided by the cryoprobes manufacturers [1]. However, in reality the combined effect of the multiple probes produces a final iceball which is very different from the union of all theoretical shapes. This approach also misses the effects of the surrounding warming structures that deform the iceball, such as large vessels. In the past years, simulations of the thermal propagation inside soft tissue for cryoablation have been investigated by a few groups. Most of the approaches implement an iterative algorithm to model the thermal propagation with a cube surrounding the tumor. The implementation is usually based on Pennes’ bioheat equation [10], which is a complex but accurate representation of the phenomenon.

Magalov et al. [9] and Ge et al. [2] both proposed simplified models to simulate iceball growth in gels. More recently, Ge et al. [3, 4] applied their model to human properties, including a simulation of the influence of vascular systems during freezing. Similar models were presented with application to various organs, e.g., the prostate [6] or the lungs [7, 8]. In most of these models, the accuracy and realism of the model come at the cost of high computation times. In 2011, Rieder et al. [11] described a GPU-based implementation of heat propagation for radio frequency ablation procedures, which is another popular type of thermal ablation using extreme heat. Their model, designed for a single probe, is based on weighted distance fields, and has been simplified to allow real-time computation. More recently, we proposed a formulation of the propagation of cold in the tissue, and an efficient implementation on GPU within a cube of interest [5].

Fig. 1.
figure 1

Examples of residual freezing along the cryoprobe’s body on two patients images, with respectively one (left) and two (right) visible cryoprobes. The dotted white lines represent the theoretical ellipsoidal shapes of the iceballs, while the red arrows highlight the residual freezing (Color figure online).

Full size image

All of the above cited approaches consider that the localization of the production of cold is limited to the freezing shaft, and that the body of the needle, supposed to be insulated, is completely passive. For all those models, the source of cold is then limited to the tips of the cryobrobes. However the well-known phenomenon, sometimes called “comet-tail effect” by radiologists, and that is very well visible on the intra-operative images, contradicts this theory. The actual dark shape of the iceball, as illustrated on Fig. 1, is not perfectly ellipsoidal but rather drop-shaped, elongated along the cryoprobe’s body. This effect is mostly due to a residual freezing along the body of the cryoprobe that is imperfectly insulated. An omission of this effect can lead to an underestimation of the predicted iceball during the planning and potential damages to healthy tissue.

In this paper, we present a method improving our previous approach to account for the residual freezing. In Sect. 2, we present our method to model the decreasing freezing along the cryoprobe, and how it is included in the simulation. In Sect. 3, we explain the experiment we conducted to validate our approach and detail the results before discussing and concluding.

2 Materials and Methods

In our previous work [5], we proposed an approach to simulate the thermal propagation within a cubic grid of voxels of 10 cm centered on the tumor, where some voxels were labeled as freezing or cooling structures with a fixed and constant temperature assigned to each label. The sources of cold were limited to the active tips of the cryoprobes.

2.1 Modelling of the Cryoprobe’s Body Temperature

To model the residual freezing along the probes, we propose an enhanced variant. We label voxels as being part of one of the three categories: (1) freezing elements (including cryoprobe’s body), (2) cooling structures and (3) normal tissue. But instead of associating a temperature to a whole set of voxels sharing the same label, we associate to each individual voxel of categories 1 or 2 its own operating temperature. This approach allows to have different temperatures, constant over time, within the same structure.

Fig. 2.
figure 2

Part of the cubic grid in which the simulations are computed. The voxels are labeled “1” if they are located along the cryoprobe. The blue color represent the freezing temperature \(T_f\), red represents the body temperature \(T_b\), and the color gradient represents the residual temperature gradient along the cryoprobe’s body. Voxels of the freezing shaft are represented with a yellow border. (Color figure online)

Full size image

The voxels of the freezing shaft are labeled “1” and homogeneously associated with the freezing temperature \(T_f\) of the corresponding cryoprobe model given by the manufacturer. The voxels located inside the vessels are labeled “2” and homogeneously associated with the body temperature \(T_b\).

The voxels located in the body of the cryoprobe are labeled “1”, and assigned an increasing temperature along its shape towards the entry point, as shown on Fig. 2.

Let us denote L the length of the part of the cryoprobe’s body that is emitting residual freezing, \(p_f\), \(\delta _f\), and \(T_f\) are respectively the center, the length, and the temperature of the freezing shaft, \(p_b\) the point along the cryoprobe’s body where the residual freezing stops and the cryoprobe’s body starts being at body temperature \(T_b\), and \(\overrightarrow{v}\) is the outwards direction vector of the cryoprobe. The notations are illustrated on Fig. 3. Temperature T at an intermediary point p along the cryoprobe is given by Eq. 1:

$$\begin{aligned} T = \lambda .T_f + (1-\lambda ).T_b \end{aligned}$$
(1)

where \(\lambda = l/L\), with \(l = \left\| p_b - p\right\| \), and \(p_b=p_f+(L+\frac{\delta _f}{2}).\overrightarrow{v}\).

The voxels labeled “0” as part of normal tissue are assigned \(T_b\) as an initial temperature, and considered as variable. The temperature of all other voxels is considered as constant.

Once built, the tridimensional matrix of labels and temperatures is passed to the thermal propagation simulation algorithm to compute the variation of temperatures over time. The advantage of assigning the constant temperatures as an initial step in a grid is that increasing the sources of cold does not have any significant influence over the computation time.

After this initial step, we use an iterative simulation algorithm that can be found in the literature, based on a discrete approximation of Pennes heat equation used in a finite space:

$$\begin{aligned} T^{new}_{i,j,k} = T_{i,j,k} + \frac{\varDelta t. \beta }{C_{i,j,k} (\varDelta x)^3} . H_{i,j,k} + C_b . \omega _b(T_{i,j,k}) . (T_b - T_{i,j,k}) + Q_m \end{aligned}$$
(2)

where \(T^{new}_{i,j,k}\) denotes the new temperature to be computed at voxel (i, j, k) after a time step of \(\varDelta t\), \(T_{i,j,k}\) and \(C_{i,j,k}\) are respectively the current temperature and the volumetric heat capacity at the same voxel, \(\beta \) is a relaxation factor set to 1.95, \(H_{i,j,k}\) is the heat flow computed from the six neighbouring cells spaced of \(\varDelta x\) from the current voxel, \(C_b\) denotes the heat capacity of blood, \(\omega _b\) is the blood perfusion rate, and \(Q_m\) is the metabolic heat rate of tissue. We refer the reader to [5] for a detailed description of the algorithm and the parameters used.

2.2 Placement of the Cube of Interest

It is important to ensure that the simulated residual freezing is fully included within the cube in which the simulation is performed. If the cube of interest is centered on the center of gravity of the cryoprobe’s tips, the shape of the residual freezing might be too close to the borders of the cube where boundary conditions apply, which would alter the simulation, then reducing and deforming the final shape.

Fig. 3.
figure 3

The drop-shaped iceball (gray) is elongated along the cryoprobe’s body. Point \(p_f\) is the center of the active freezing shaft (yellow). A residual freezing is observed over a length l of the body of the cryoprobe, up to point \(p_b\). Intermediary point \(p_i\) produces residual freezing at a temperature \(T_i\). (Color figure online)

Full size image

One possibility is simply to enlarge the cube enough to ensure a good coverage of the whole final iceball. However, this approach is very time consuming as an increase of \(10\%\) of the size of the cube implies an increase of \(33\%\) of the computation time. If the computation time, and thus the reduction of the size of the cube, is an important issue, then a good centering of the final iceball including the residual freezing will be essential to ensure a fast and accurate computation and avoid simulations in irrelevant parts of the volume.

In order to avoid this effect, it is recommended to shift the center of the cube along the vector \(\overrightarrow{V}\) averaging the direction vectors of all cryoprobes. For N cryoprobes, if \(\overrightarrow{v_n}\) is the direction vector of the \(n^{th}\) cryoprobe, \(\overrightarrow{V}=\frac{\sum _{n=1}^{N} \overrightarrow{v_n}}{N}\).

The amount of the shift will depend on several factors, among which the size of the cube of interest, the chosen value of L, the cryoprobe model (on which depends the size of the iceball), and the global orientation of the cryoprobes (towards a corner or a side of the cube). In this work, we considered a cube of 8 cm to be in most cases a good compromise.

3 Experiment and Results

We experimented our simulation on 5 datasets of patients who underwent renal cryoablation interventions. For all interventions, intraoperative T1 MR images were acquired during the maximal freezing session, showing quite clearly the hypodense signal corresponding to the iceball generated from the cryoprobes. The interventions were performed using 3 to 5 cryoprobes of two models: IceRod and IceEdge from Galil Medical, interacting to produce a large iceball covering the whole tumor. Each patient case had only one single model of probes. The dimensions of the probes are respectively: IceEdge 2.4 mm (10G) thick, with an active freezing part of 28 mm starting at 5.2 mm from the tip, and IceRod 1.5 mm (17G) thick, with an active freezing part of 31 mm starting at 4.2 mm from the tip.

Fig. 4.
figure 4

Result of the simulation on cases 3 (left) and 4 (right): simulated iceball shapes with residual freezing \(\mathcal {I}_{RF}\) (blue) and without \(\mathcal {I}_{ELL}\) (green), and segmented shape \(\mathcal {I}_{SEG}\) (red). (Color figure online)

Full size image

For each dataset, the iceball including the residual freezing (\(\mathcal {I}_{SEG}\)) was interactively segmented and the entry and tip points of each cryoprobe were manually selected under the supervision of an expert radiologist. The simulation was computed twice: once using our approach considering the residual freezing (RF), and once using the approach from the literature [5] without residual freezing (ELL), producing the respective resulting shapes \(\mathcal {I}_{RF}\) and \(\mathcal {I}_{ELL}\) corresponding to the \(0^{\circ }\)C isotherm surfaces. Simulation with no residual freezing was achieved by labelling “0” the cryoprobes bodies. All the simulations were run on a desktop computer equipped with a core-i5 3.20 GHz CPU with 32 Gb RAM and a GeForce GTX-760 GPU.

In all cases a cube of 8 cm was used to allow space around the iceball and avoid the direct influence of the boundary condition. Length L has been chosen empirically according to the cryoprobes types and their respective size: 45 mm for IceRod, and 60 mm for IceEdge. A shift of \(15\%\) of L has been used in all cases.

Each experiment simulated a full standard procedure with alternating freezing/thawing cycles: 10 min of active freezing, followed by 9 min of passive thawing and 1 min of active thawing, and again 10 min of active freezing. We used the parameters provided by the manufacturer: the freezing temperature at the cryoprobe’s tip was set to \(-138.0^{\circ }\)C for IceEdge, and to \(-119.40^{\circ }\)C for IceRod. Results are illustrated on Fig. 4.

To evaluate the accuracy of our simulation, \(\mathcal {I}_{RF}\) and \(\mathcal {I}_{ELL}\) were compared to \(\mathcal {I}_{SEG}\) by computing the Hausdorff distance and the Dice coefficient. The results, along with computation times, are presented on Table 1.

Table 1. Dice coefficient and Hausdorff distance (HD, in mm) measurements between the computed and the ground truth iceball segmentation shapes, and computation times, for simulations with and without residual freezing.
Full size table

It can be observed that the simulated shape \(\mathcal {I}_{RF}\) fits the segmented shape \(\mathcal {I}_{SEG}\) better than \(\mathcal {I}_{ELL}\). The average Dice coefficient is 0.821, and the average Hausdorff distance is 9.340 mm, improving the results of the simulation without residual freezing. As expected, the computation time is not significantly different when considering a simulation with or without residual freezing. The computation time mostly depends on the size of the cube, i.e. the number of voxels to explore during the simulation.

This first attempt of simulating the residual freezing has been designed assuming that the increase of temperature along the cryoprobe’s body was linear. It would be interesting to perform experimental measurements to validate this theory, and get more accurate parameters to refine our model. Non linearity of the temperature distribution could be achieved by replacing Eq. 1 by a more complex model.

In our experiments, we chose to use a cube size of 8 cm, as a trade-off between a full inclusion of the simulated iceball within the cube with enough distance to the boundaries, and a reduced computation time. For case 5, a smaller cube of 6 cm could have been used, as the number and model of cryoprobes are smaller. However, it could be interesting to investigate the automatic suggestion of an optimal cube size, based on the number of cryoprobes, their orientation in the image, and their model.

Previous works in the literature obtained simulation results that were usually compared against theoretical shapes of iceballs (pure ellipsoids) or segmented iceballs in which the residual freezing part was usually not included. Yet, modelling the residual freezing is clinically relevant and important, as an underestimation of the predicted iceball shape may result in potential damage to healthy cells, causing pathologies or pain. With the objective in mind to include such a simulation within a preoperative treatment planning framework, with an automatic computation of optimal placements for the cryoprobes, an accurate and complete prediction of the iceball is essential.

4 Conclusion

In this paper we presented a first approach to model the residual freezing along the cryoprobes bodies during renal cryoablation. This phenomenon was usually ignored in previously published works, although it may result in an underestimation of the predicted frozen area and to potential damage to healthy tissue, or cause pain. Our approach consists in modelling a temperature linearly increasing along the cryoprobe’s body but constant over time, and use it as a source of cold during the simulation. We experimented our model on 5 retrospective patient datasets, and showed that our model fits the segmented iceball with a reasonable accuracy. We determined that a size of 8 cm for the simulation cube, with an appropriate shift, is enough to simulate the whole 30 mn process and produce the shape of the resulting iceball in less than 20 s. Further experiments should be performed in future works to confirm these first results.