CC BY
114
ADVANCED RESEARCHES
Calculation of Effective Freezing Time
in Lung Cancer Cryosurgery Based on Godunov Simulation
DOI: 10.17691/stm2016.8.1.07 Received December 24, 2015
EG. Kotova, MD, PhD, Researcher1- 2;
V.I. Kochenov, MD, DSc, Senior Researcher, Department of Operative Surgery and Topographic Anatomy1; Director2;
S.N. Tsybusov, MD, DSc, Professor, Head of the Department of Operative Surgery and Topographic Anatomy, Vice-Rector for Academic Affairs1;
D.Y. Madai, MD, DSc, Professor, Head of the Department of Maxillofacial Surgery and Surgical Dentistry, Faculty of Dentistry and Medical Technology3;
AV. Gurin, Surgeon, Oncologist, Clinical and Diagnostic Department4
1Nizhny Novgorod State Medical Academy, 10/1 Minin and Pozharsky Square, Nizhny Novgorod, 603005, Russian Federation;
Scientific Clinical Center Medical Cryology “onKolor”, 4a Beketova St., Nizhny Novgorod, 603057,
Russian Federation;
3Saint Petersburg State University, 7-9 Universitetskaya naberezhnaya, Saint Petersburg, 199034,
Russian Federation;
4Petrov Research Institute of Oncology, Ministry of Health of Russia, 68 Leningradskaya St., Pesochnyy Village, Saint Petersburg, 197758, Russian Federation
There have been presented the results of lung cancer cryosurgery simulation using numerical solutions of enthalpy equation according to Godunov method. For the cryodestruction improvement purposes we successfully calculated the effective freezing time taking into account the evolution of an ice ball covering the tumor area. Geometrical transformation parameters of an ice ball have been measured by calculating the temperature distribution and the interface position in biological tissue. Mathematical cryosurgical procedures are described by heat transfer equations in solid and liquid phases. Numerical results for one-dimensional case were verified by comparing with exact solutions. In twodimensional modeling an effective cryotherapy time, which corresponds to freezing time of all tumor parts, was calculated as the area of forming ice balls covering all tumor region. The findings enable to set the effective time of a cryosurgical procedure in lung cancer. The knowledge of temperature distribution and interface position in biological tissue offers an opportunity to a cryosurgeon to finish the procedure within a certain time period to minimize the healthy tissue damage and destroy tumor cells to the maximum. Simulation application enables to schedule cryotherapy in lung cancer more effectively and to a good quality.
Key words: cryosurgery; lung cancer; Godunov method; cryosurgery simulation.
Lob-, bilob- and pneumonectomy with ipsilateral mediastinal lymph node dissection is a conventional treatment method for stage I-II and ША-ШВ 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
For contacts: Tatiana G. Kotova, e-mail: 261126.79@mail.ru
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
48 СТМ J 2016 — vol. 8, No.1 T.G. Kotova, V.I. Kochenov, S.N. Tsybusov, D.Y. Madai, A.V. Gurin
ADVANCED RESEARCHES)
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 Q 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 QS and QL 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 x=(x,y)eR2 and time t. In the frozen area (solid phase) the heat transfer equation can be expressed as follows:
p c 2TM - k (0"T
rscs
dx2
x,t) ^ d2T | x
’*)
Эу2
XEQs, (1)
where ps, cS and kS 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.
Plcl
dT (x,t) ~dT~
( Д2
=kL
d2T (
dx2
x,t ) + d2T i x,
,t)'
Эу2
(2)
+ ШЬPbcb [Tb-T (X,t)] + Xe^L-
where pL, cL and kL is density, specific density and heat transfer of unfrozen tissue, respectively; wb, pb, cb and Tb is perfusion, density, specific heat capacity and blood temperature, respectively, Qm is m eta bolic h eat evol uti on. Analytic studies of bioheat transfer sire given in Hterature
[12, 13]. In the present study it is accepted that pS=pL=p, therefore, during the freezing process the volume does not expand.
The condition of interface position corresponds to Stefan condition:
PLvn =- kL
dT (x+ ,t) Эп
+ ks
dT (x,t) Эп
хеГ,
(3)
where L is latent heat, vn 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:
T(x,t)=Tm, хеГ, (4)
where Tm 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:
E (x,t)
pcs (T (x,t)~ Tm ), T ( x,t )< Tm ;
PCL (T (x,t)~ Tm ) + PL, T (x,t)> Tm ,
where T(x,t)<Tm and T(x,t)>Tm are temperatures of solid and liquid phase, respectively.
Suppose 0^x^l1, 0^y^l2 is two-dimensional region of biological tissue, in this case, l1=l2=0.4 m. The region [0, l1], [0, l2] is divided into subintervals M1 and M2, respectively. Thus, we obtain control volumes M1 and M2. The internal area Vu=[xi_1l2, x,+1/2]x[y_1/2, y+1/2] is defined as control volume, where x-1/2 is the node between x-1 and x. The conservation of energy in each control volume Vu can be expressed as
J [E (x,t+At )-E (x,t)] dA=J J -qn dSdt,
(6)
where E(x,t) is the enthalpy per unit area, -qn is heat flux into the area ^across its boundary dVj n being the outgoing unit normal to dVu.
An explicit scheme in two-dimensional region (6) based on Godunov method will be represented as
ЕП + l _ ГЛ .
Eij - 4,J + ДХ
[q/-l/2,J qi+l/2,J ] +
jy [air-V2 - air+|/2 ] + AHPbcb [ Tb -Ti ] +AQm 1
(7)
Calculation of Effective Freezing Time in Lung Cancer Cryosurgery Based on Godunov Simulation СТМ J 2016 — vol. 8, No.1 49
ADVANCED RESEARCHES
where
qi -1/2,j ~
_ - tu
Ax
, Ri-1/2,j - 2
4-1/2,j p ’ ■ i-1/2,j 9
rV-1/2,j 2
q _ Ti,j -1 ~ Ti,j R _АУ
4i,j-1/2 _ p ’ -1/2 _ 9
Ri,j-1/2 2
1 1
K/-u + ku
11
V ki,j-1 ki,j J
Temperature distribution in tumor cells and healthy tissues is derived from the equation (5), where enthalpy value in each control volume Vi; 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 TL>Tm. As x=0 temperature is maintained at the level TS<Tm, freezing starts in the direction from the left to right margin of the region, if 0<x<r(t), and Г<x«» is the regions of solid and liquid phases, respectively. Exact solutions for this problem are the following:
T (xt ) =
erf
T+(Tm - T )-
ЩТ1
erf (A)
, 0<x<Г (t);
erfc
T -( T - Tm
л/2оА
(8)
erfc (KyjaJaL
?x>Г(t)
and
Г (t)= 2\JaSt, t>0.
(9)
The term erf denotes error function, and parameter l is obtained by solving a transcendental equation:
SL
StL
where k
exp (A2) erf (A) exp (v2 Л2) erf (vA)
Cs T -Ts).
= W^, (10)
as=-K-; aL=-K-; Sts = CsTL s PCs L POl s L
X
X
StL=Cl T Tm ^; v=-JaJaL.
Godunov method for one-dimensional case is represented 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 Г is approximated:
r"=Xm-V2 +(1-K,)*. (11)
Here m is an index, wherein control volume Vm contains interface, and lm is liquid phase, which can be expressed as:
К
0 , E- <0 (s o I i d f r actio n) ; En
-j-, 0<En<pL (interface); 1, ЕП >pL (liquidfraction).
(12)
Figures 1 and 2 show the temperature distribution and interface position calculated by Godunov method, and their comparison with exact solutions [See (7)-(12)].
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 Ax=0.1/320 equal to 1.76 and 0.013%, respectively.
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°C; TS=-196°C; length /=0.1 m; Ax=0.1/320
Figure 2. Interface position r(t) for one-dimensional case considering the time; calculated using physical properties (See the Table); TL=37°C; TS=-196°C; length /=0.1 m; Ax=0.1/320
///////////////////////////////////////////////////////////////^^^^
50 СТМ J 2016 - vol. 8, No.1 T.G. Kotova, V.I. Kochenov, S.N. Tsybusov, D.Y. Madai, A.V. Gurin
ADVANCED RESEARCHES
Physical properties of tissues
Signs Parameter Value Unit of measurement
Cs Specific heat capacity of frozen lung tumor tissue 1.23 kJ/kg/°C
Cl Specific heat capacity of unfrozen lung tumor tissue 4.2 kJ/kg/°C
Cb Blood specific heat capacity 3.64 kJ/kg/°C
ks Heat transfer of frozen lung tumor tissue 2.2510-3 kJ/kg/°C
kL Heat transfer of unfrozen lung tumor tissue 0.5510-3 kJ/kg/°C
Tm Melting temperature 0 °C
L Latent heat 333 kJ/kg
P Lung tumor tissue density 1,000 kg/m3
Pb Blood density 1,000 kg/m3
Qm Metabolic heat generation in a tumor 42 kJ/s/m3
Wb Blood perfusion in a tumor 0.002 ml/s/ml
Tb Blood temperature 37 °С
To Initial temperature 37 °С
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).
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).
The initial area of tumor region and the left lung was 0.001899 and 0.035359 m2, respectively. By the time f=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).
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
Figure 3. Lung cancer scheme, with an area exposed to cryosurgery marked
Calculation of Effective Freezing Time in Lung Cancer Cryosurgery Based on Godunov Simulation СТМ J 2016 — vol. 8, No.1 51
ADVANCED RESEARCHES
0.Z5
u.n s
tm 112.92 s
t- 211.72
CT)
CD
X
0.2
0 15
0.Z5
CT)
CD
X
0-2
| FH - D.№(l№ m^.lSSi I FC - 0 000*S6 rii* (26.1144
I - 310 52 s
I FH - HOW7% m7{2.2$4)
I FC - Q.Q0l?ro m* [-Sd-Ei-’^h
(MS'—
0.3
0 35
I FH - OW10S7 П1 I FC - 0.0016£? ip3 (W.OiTt)
I FH - Q OOHje irt1 |рС-й№19Э0л^(Ю0.№|
0.4 0.3
0.35
0.4 0.3
0.35 0.4
Width (m)
s
Figure 4. The evolution of ice balls during freezing, with percentage ratio of the areas of frozen healthy tissue (FH) and the area of frozen tumor cells (FC). The sequence of the images is shown from the upper left to the lower right
Figure 5. The position of six chosen points with an ice ball on the scheme of the lung cancer region under study
Figure 6. Temperature changes at the chosen six points during freezing
«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 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.
Thus, the present study succeeded in calculating effective freezing time in lung cancer cryosurgery taking into consideration the evolution of an ice ball covering
52 СТМ j 2016 — Vol. 8, No.1 T.G. Kotova, V.I. Kochenov, S.N. Tsybusov, D.Y. Madai, А .V. Gurin
ADVANCED RESEARCHES
the 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.
References
1. Spravochnik po onkologii [Handbook of oncology].
Pod. red. Moiseenko V.M. [Moiseenko V.M. (editor)].
Saint Petersburg; 2008.
2. Kumar S., Katiyar V.K. Numerical study on phase change heat transfer during combined hyperthermia and cryosurgical treatment of lung cancer. Int J of Appl Math and Mech 2007, 3(3): 1-17.
3. Wan R., Liu Z., Muldrew K., Rewcastle J. A finite element model for ice ball evolution in a multi-probe cryosurgery. Comput Methods Biomech Biomed Engin 2003; 6(3): 197-208, http:// dx.doi.org/10.1080/1025584031000151185.
4. Yang B., Wan R.G., Muldrew K.B., Donnelly B.J. A finite element model for cryosurgery with coupled phase change and thermal stress aspects. Finite Elem Anal Des 2008; 44(5): 288297, http://dx.doi.org/10.1016Zj.finel.2007.11.014.
5. Rossi M.R., Tanaka D., Shimada K., Rabin Y. An efficient numerical technique for bioheat simulations and its application to computerized cryosurgery planning. Comput Methods Programs Biomed 2007; 85(1): 41-50, http://dx.doi.org/10.1016/ j.cmpb.2006.09.014.
6. Rossi M.R., Tanaka D., Shimada K., Rabin Y. Computerized planning of cryosurgery using bubble packing: an experimental validation on a phantom material. Int J Heat and Mass Transf 2008; 51(23-24): 5671-5678, http://dx.doi. org/10.1016/j.ijheatmasstransfer.2008.04.045.
7. Shi J., Chen Z., Shi M. Simulation of heat transfer of biological tissue during cryosurgery based on vascular trees. Applied Thermal Engineering 2009; 29(8-9): 1792-1798, http:// dx.doi.org/10.1016/j.applthermaleng.2008.08.014.
8. Zhao G., Zhang H.-F., Guo X.-J., Luo D.-W., Gao D.-Y.
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
Effect of blood flow and metabolism on multidimensional heat transfer during cryosurgery. Med Eng Phys; 2007; 29(2): 205215, http://dx.doi.org/10.1016/j.medengphy.2006.03.005.
9. Shih T.-C., Yuan P., Lin W.-L., Kou H.-S. Analytical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface. Med Eng Phys 2007; 29(9): 946-953, http://dx.doi.org/10.1016/j.medengphy.2006.10.008.
10. Chua K.J., Chou S.K., Ho J.C. An analytical study on the thermal effects of cryosurgery on selective cell destruction. J Biomech 2007; 40(1): 100-116, http://dx.doi.org/10.1016/ j.jbiomech.2005.11.005.
11. Deng Z.-S., Liu J. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. J Biomech Eng 2002; 124(6): 638-649, http://dx.doi.org/10.1115/1.1516810.
12. Chua K.J., Chou S.K., Ho J.C. An analytical study on the thermal effects of cryosurgery on selective cell destruction. J Biomech 2007; 40(1): 100-116, http://dx.doi.org/10.1016/ j.jbiomech.2005.11.005.
13. Voller V.R., Shadabi L. Enthalpy methods for tracking a phase change boundary in two dimensions. International Communications in Hegt&nd Mass Transfer 1984; 11(3): 239249, http://dx.doi.o^50i016/0735-1933(84)90040-x.
Calculation of Effective Freezing Time in Lung Cancer Cryosurgery Based on Godunov Simulation СТМ J 2016 — vol. 8, No.1 53
