Calculation of Effective Freezing Time in Lung Cancer Cryosurgery Based on Godunov Simulation

Lob-, bilob- and pneumonectomy with ipsilateral mediastinal lymph node dissection is a conventional treatment method for stage I–II and IIIА–IIIВ non-small cell lung carcinoma. Currently, in operable cases, there has been studied the efficiency of both neoadjuvant and adjuvant chemotherapy. Radiotherapy is a method of choice in inoperable cases. Palliative chemotherapy is chosen in stage IV lung cancer. Small cell lung cancer treatment includes various combinations of chemo- and radiotherapy. Surgical technique is considered as a part of combined treatment in patients with stages I and II [1].

Cryosurgery is one of surgical techniques, when extremely low temperatures are used to destroy tumor tissues. In recent years, there have been developed cryosurgical techniques to treat severe cancer types, such as brain cancer, breast cancer, prostate carcinoma, renal cell carcinoma, hepatic cancer. Tumor cells are exposed to extremely low temperatures (liquid nitrogen at –196°C) using a cryoprobe. Due to low temperatures, an ice ball forms around a cryoprobe, the ice ball keeping continuously freezing tumor cells. Cryotherapy results in biological tissue separation into two parts: the solid substance and the liquid one. Frozen tumor cells are damaged if their temperature is below –30°C [2]. The procedure aims at maximum damage to malignant tumor cells with minimal damage to surrounding healthy tissue.

Difficulties in cryosurgery stem due to the necessity to minimize the healthy tissue damage that has determined active studies of a freezing process during a cryosurgical procedure using computerized simulation. Wan et al. [3] have demonstrated the process of the ice ball evolution during a cryosurgical procedure using a finite-element method for its analysis. Similar numerical method was also used in the work [4] to simulate prostate cryosurgery, the aspects of thermal stress being considered. Moreover, Rossi et al. [5, 6] developed an effective numerical method of computer-aided planning for cryosurgery. Other examples of simulation in cryosurgery can be found in the researches by Shi and Zhao [7, 8]. However, these models do not take into consideration an effective freezing time regardless the fact that a time factor is an essential aspect to attain success when carrying out cryosurgical procedures.

The aim of the investigation was to determine an effective freezing time in lung cancer cryosurgery using the simulation of interface position resulted from temperature distribution by studying the ice ball evolution.

Cryosurgical processes are simulated mathematically as the equations of heat transfer in solid and liquid phases, where the interface region between two phases is subject to Stefan condition. In particular, the equation of heat transfer in liquid phase is presented as Pennes bioheat transfer equation [9–11]. Due to the difficulty of equation solution, the equations are reformulated as energy equations (enthalpy). The advantage of applying enthalpy calculation is in fact that basic equations are unaltered regardless a phase they are applied to: liquid or solid, therefore, for their solution standard numerical schemes can be easily applied, e.g. Godunov method [12, 13].

Simulation mathematical tools

The region Ω under study is a part of tissue occupied by healthy and tumor cells, which initially are in liquid state. When tumor cells are started freezing, there is the phase change from liquid into solid.

Let Ω_S and Ω_L be the region of solid and liquid phase, respectively, and Г is clear and soft interface separating the areas of solid and liquid states. Suppose that T(x,t) is the temperature at position and time t. In the frozen area (solid phase) the heat transfer equation can be expressed as follows:

where ρ_S, c_S and k_S is density, specific density and heat transfer of frozen tissue, respectively.

However, in an unfrozen area (liquid phase) due to blood perfusion and metabolic process, the heat transfer equation can be written in the form of a bioheat transfer equation, i.e.

where ρ_L, c_L and k_L is density, specific density and heat transfer of unfrozen tissue, respectively; ω_b, ρ_b, c_b and T_b is perfusion, density, specific heat capacity and blood temperature, respectively, Q_m is metabolic heat evolution. Analytic studies of bioheat transfer are given in literature [12, 13]. In the present study it is accepted that ρ_S=ρ_L=ρ, therefore, during the freezing process the volume does not expand.

The condition of interface position corresponds to Stefan condition:

where L is latent heat, n_n and n are normal components of velocity, and an outgoing unit perpendicular to Г, respectively.

Moreover, the temperature on interface surface can be expressed as follows:

where T_m is melting temperature.

Godunov method implementation

Since interface position is unknown and should be determined for each case, numerical calculation of the equations (1)–(4) is indirect. The techniques to solve the problem will consist in the rearrangement of heat transfer equations of solid and liquid phases into energy equation (enthalpy). In enthalpy, the boundary between solid and liquid phases can be disregarded, so that a numerical scheme can be easily applied in terms of energy conservation, i.e. Godunov method.

Suppose E(x,t) denotes the enthalpy per unit area at position x and time t, the sum of sensible heat and latent heat is:

where T(x,t)<T_m and T(x,t)>T_m are temperatures of solid and liquid phase, respectively.

Suppose 0≤x≤l₁, 0≤y≤l₂ is two-dimensional region of biological tissue, in this case, l₁=l₂=0.4m. The region [0, l₁], [0, l₂] is divided into subintervals M₁ and M₂, respectively. Thus, we obtain control volumes M₁ and M₂. The internal area V_i,j=[x_i–1/2, x_i+1/2]×[y_i–1/2, y_i+1/2] is defined as control volume, where x_i–1/2 is the node between x_i–1 and x_i. The conservation of energy in each control volume V_i,j can be expressed as

where E(x,t) is the enthalpy per unit area, –qn is heat flux into the area V_i,j across its boundary ∂V_i,j, n being the outgoing unit normal to ∂V_i,j.

An explicit scheme in two-dimensional region (6) based on Godunov method will be represented as

where

Temperature distribution in tumor cells and healthy tissues is derived from the equation (5), where enthalpy value in each control volume V_i,j is calculated using the equation (7).

The model verification for one-dimensional case

Exact solutions of the equations (1)–(4) are available for semi-infinite one-dimension cases, though without taking into account blood perfusion and metabolic heat generation. The structure of one-dimensional problems in cryosurgery can be described as follows. Suppose that 0<x<∞ is semi-infinite domain, which initially is in the liquid phase at temperature T_L>T_m. As х=0 temperature is maintained at the level T_S<T_m, freezing starts in the direction from the left to right margin of the region, if 0<x<Г(t), and Г<x<∞ is the regions of solid and liquid phases, respectively. Exact solutions for this problem are the following:

and

The term erf denotes error function, and parameter λ is obtained by solving a transcendental equation:

where

Godunov method for one-dimensional case is repre­sented by the equation (7), but without two last members in the right side of the equation. The temperature was calculated numerically using the equation (5), with interface surface position at a time point tⁿ is approximated:

Here m is an index, wherein control volume V_m contains interface, and λ_m is liquid phase, which can be expressed as:

Figures 1 and 2 show the temperature distribution and interface position calculated by Godunov method, and their comparison with exact solutions [See (7)–(12)].

---

Figure 1. Temperature distribution T(x,t) for one-dimensional case at time t=522.24 s; calculated using physical properties (See the Table); TL=37°С; TS=-196°С; length l=0.1 m; Δx=0.1/320

---

Figure 2. Interface position Г(t) for one-dimensional case considering the time; calculated using physical properties (See the Table); TL=37°С; TS=-196°С; length l=0.1 m; Δx=0.1/320

---

Physical properties of tissues

---

The figures clearly demonstrate that numerical Godunov scheme agrees completely with an exact solution. A mean error of temperature distribution and interface position (with the assistance of physical properties, See the Table) at Δx=0.1/320 equal to 1.76 and 0.013%, respectively.

Computerized lung cancer cryosurgery simulation

Simulation models are designed to help a surgeon to determine the duration of lung cancer cryosurgery procedure. As already noted, time factor in this case is of great importance, since it enables to reduce the damage risk of healthy tissues preventing them from freezing. We are going to demonstrate it with an example of the left lung human cancer of randomly chosen geometrical forms and location, and show how an effective freezing time is calculated (Figure 3).

---

Figure 3. Lung cancer scheme, with an area exposed to cryosurgery marked

---

Extremely low temperature (–196°С) is given through a cryoprobe to tumor cells, which are initially in liquid phase at 37°С. When the tissue temperature is falling, an ice ball forms around the cryoprobe, the ice starting extending to the outside of the cryoprobe, into tumor cells surrounding them. At certain time an ice ball covers all target areas. Effective freezing time is considered to be the time, within which all target areas containing tumor cells are frozen.

To determine effective freezing time, it is necessary to calculate the area of frozen tumor cells and healthy tissue considering the ice ball evolution. To simulate changing geometrical forms of an ice ball covering the tumor area we used the physical properties given in the table below.

Figure 4 shows the process of ice ball evolution with the percentage ratio of the area of frozen tumor cells (FC) and frozen healthy tissue (FH) on the region under study (See Figure 1).

---

Figure 4. The evolution of ice balls during freezing, with percentage ratio of the areas of frozen healthy tissue (FН) and the area of frozen tumor cells (FC). The sequence of the images is shown from the upper left to the lower right

---

The initial area of tumor region and the left lung was 0.001899 and 0.035359 m², respectively. By the time t=14.11 the ice ball occupied 26.11% of the area of tumor cells and 0.18% of the healthy tissue area. Tumor cells were covered completely by the time t=522.24 s, the coverage area of frozen healthy tissue being 4.06%. Therefore, cryosurgical procedure is to be stopped at the time t=522.24 s to prevent healthy tissue damage growth.

For temperature recording, on the scheme we have chosen six points inside and outside the lung cancer region (Figure 5).

---

Figure 5. The position of six chosen points with an ice ball on the scheme of the lung cancer region under study

---

It can be noticed that the point 1 situated near to the cryoprobe was frozen within less than 7 s, while points 2 and 3 needed 84 s to become frozen. Freezing was completed 522.24 s after the process initiation; the temperature at points 1, 2 and 3 being –140, –103 and –67°С, respectively. Thus, tumor cells around these three points were damaged. At point 4 unwanted freezing of healthy tissue occurred by the time of 331 s, and by the end of the exposure the temperature was ≈19°С. Moreover, it should be noted that 522.24 s after the procedure started there was no freezing at points 5 and 6. It corresponds to the cryosurgery purposes: these points are beyond the tumor region and therefore, they were not to be exposed to freezing.

Thus, the first three points are damaged, while the last three points have normal temperature of a healthy human. Figure 6 shows the temperature changes at these six points.

---

Figure 6. Temperature changes at the chosen six points during freezing

---

Figure 7 shows the temperature distribution during the cryosurgery procedure for several time points. At t=522.24 s, nearly all tumor parts had the temperature below –50°С, as a result, the cancer cells inside that area were damaged. The figure also reveals the interface position, which describes the geometry of the ice ball.

---

Figure 7. Temperature distribution and interface position during a cryosurgical procedure. The sequence of the images is shown from the upper left to the lower right

---

Thus, the present study succeeded in calculating effective freezing time in lung cancer cryosurgery taking into consideration the evolution of an ice ball coveringthe tumor area. Geometrical parameters of ice ball evolution were obtained by studying temperature distribution and interface position in biological tissue. Effective time for cryosurgical procedure was found to be 8 min 42 s. This means that in similar cases of lung cancer a cryosurgeon is to finish the procedure at this time to prevent healthy tissue from damaging. The knowledge of temperature distribution and interface position in biological tissue will enable to minimize the healthy tissue damage and destroy tumor to the maximum.

Conclusion. The calculation of effective freezing time and cryosurgery simulation based on the calculation enable to plan cryosurgery in lung cancer more effectively and to a good quality.

Study Funding and Conflicts of Interest. The study was not funded by any sources, and the authors have no conflicts of interest related to the present study.