Fast inverse prediction of the freezing front in cryosurgery

Author’s Accepted Manuscript Fast inverse prediction of the freezing front in cryosurgery Mohamed Hafid, Marcel Lacroix

www.elsevier.com/locate/jtherbio

PII: DOI: Reference:

S0306-4565(17)30132-8 http://dx.doi.org/10.1016/j.jtherbio.2017.05.008 TB1940

To appear in: Journal of Thermal Biology Received date: 6 April 2017 Revised date: 29 May 2017 Accepted date: 30 May 2017 Cite this article as: Mohamed Hafid and Marcel Lacroix, Fast inverse prediction of the freezing front in cryosurgery, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2017.05.008 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Fast inverse prediction of the freezing front in cryosurgery

Mohamed Hafid*, Marcel Lacroix1 Département de génie mécanique, 2500 boulevard de l’Université, Sherbrooke (Québec), Canada J1K 2R1.

[email protected] [email protected]

*Corresponding author. Address: Université de Sherbrooke, Département de génie mécanique, 2500 boulevard de l’Université, Sherbrooke (Québec), Canada J1K 2R1. Local : C1-4012, Tel.: +1 873 888 7768

Abstract Cryosurgery has become a well-established technique for the ablation of undesirable tissues such as tumors and cancers. The motivation for this study is to improve the efficacy and safety of this technique. This study presents an inverse heat transfer method for monitoring the motion of the freezing front from a cryoprobe. With the help of a thermocouple inserted into the layer of diseased tissue, the inverse heat transfer method estimates simultaneously the blood perfusion rate and the thermal conductivities of both frozen and unfrozen tissues. This information is then fed to the Pennes bioheat equation that: (1) calculates the time-varying temperature distribution inside the layer of tissue and (2) predicts the motion of the freezing front. The effect of the most influential parameters on the inverse predictions is investigated. These parameters are (1) the initial guesses for the unknown Levenberg-Marquardt polynomial parameters of the thermo-physical properties; (2) the temperature of the 1

Université de Sherbrooke, Département de génie mécanique, 2500 Blvd. de l’Université,

Sherbrooke (Québec), Canada J1K 2R1. Local : C1-4071, Tel. : (819) 821-8000 poste 62145

1

cryoprobe; (3) the heat transfer coefficient of the impinging jet of liquid nitrogen; and (4) the noise on the temperature data recorded by the thermocouple probe. Results show that the proposed inverse method is a promising alternative to ultrasound and Magnetic Resonance Imaging (MRI) for monitoring the motion of the freezing front during cryosurgery. For all the cryogenic scenarios simulated, the predictions of the inverse model remain accurate and stable.

Keywords: Inverse problem; Bioheat transfer; Cryosurgery; Phase change interface; Levenberg-Marquardt Method; Broyden method; Pennes equation; Enthalpy method.

NOMENCLATURE

dt fs

specific heat [J/kg°C] time step [s] solid fraction

h

heat transfer coefficient [W/m2 °C]

I

total number of measurements Jacobian matrix thermal conductivity [W/m °C] latent heat of freezing [J/kg] thickness of tissue [m]

c

J

k L Ltissue

N P Qm RRMSE

number of unknown parameters vector of unknown parameter metabolic heat generation[W/m3]

Error

relative root-mean-square error [%] estimation error [%]

s t 

freezing front position [m]

t

time [s]

Tˆ Tb

estimated temperature [°C] arterial blood temperature [°C]

Tc x

Y wb

cryogenic probe temperature [°C] Cartesian spatial coordinate [m] measured temperature [°C] blood perfusion rate [ml/s ml]

Greek symbols 2

 

 



 Pˆ

thermal diffusivity [m2/s] small number damping parameter density of tissue [kg/m3] standard deviation of the measurement error standard deviation of the

j

estimated parameter

  H  k  

sum of squares norm small number enthalpy [J/m3] difference diagonal matrix Neumann’s parameter random number

Subscripts  b exact f

ambient blood exact solution frozen region

liq max

liquidus maximum

u

unfrozen region

P

parameter

sol T

solidus temperature

Superscripts

k

T

iteration number transposed matrix estimated parameter vector matrix

3

1. Introduction Cryosurgery is the use of extreme cold in surgery to destroy abnormal or diseased tissue. Warts, moles, skin tags and small skin cancers are candidates for cryosurgical treatment. Internal disorders, including cancers, are also treated with cryosurgery. Indeed, cryosurgery is a minimally invasive procedure. It is often preferred to more traditional kinds of surgery because it minimizes pain and scarring, and its costs are relatively low (Malawer et al. 1999; Korpan, 2012). It works by taking advantage of the destructive effect of freezing temperatures on cells. When their temperature sinks beyond a certain level, ice crystals begin forming inside the cells. Due to their lower density, the ice crystals tear apart the cells. They also harm the blood vessels that supply the malignant tissue. A common method of freezing lesions is using liquid nitrogen as the cooling solution. This very cold liquid, its temperature is -196 0C, may be sprayed on the diseased tissue or circulated through a cryoprobe (Ravikumar et al. 1991; Gage and Baust, 2007). The cryoprobe is a tube placed in contact with the tumor. Ultrasound or magnetic resonance imaging (MRI) is used to guide the cryoprobe and monitor the freezing of the cells (Pfleiderer et al. 2002; Huber et al. 2001). This helps in limiting damage to adjacent healthy tissues and, more of a concern, avoiding damage to nerve tissue. But ultrasound and MRI monitoring techniques are not always readily available. The alternative then is to resort to the mathematical modeling and the numerical simulation of the heat transfer and freezing processes that take place inside the biological tissues (Jaeger and Carin, 2002; Zhao and Chua, 2012; Khanday and Hussain, 2015; Carin and Jaeger, 2001). This computational monitoring technique is convenient and relatively inexpensive. Its predictions are reliable provided, however, that the thermal properties of the tissues, such as the blood perfusion and the thermal conductivities, are well defined. Unfortunately, the thermal properties of the biological tissues are often poorly known. Moreover, they change as the tissues freeze up (Cooper and Trezek, 1971; Lee and Bastacky, 1995). It is possible to determine experimentally the thermal properties of biological tissues (Delhomme et al. 1992; Burger and Breukelen, 2013). But the measurement techniques are not without complications and they require sophisticated equipment that can hardly be employed during cryosurgery.

4

A promising approach for handling this problem rests on inverse heat transfer techniques (Ozisik and Orlande, 2000; Gosselin et al. 2009). Investigations based on inverse methods have been carried out in the past for estimating the thermal parameters of living tissue (Lesnic and Ivanchov, 2015; Huntul et al. 2017; Cao and Lesnic, 2016; Hazanee and Lesnic, 2014). Trucu et al. (2008, 2009, 2010a, 2010b, 2011), Grabski et al. (2016) and Lesnic (2009) estimated the blood perfusion coefficient for the one-dimension Pennes equation using the Tikhonov regularization method. Lee et al. (2013) determined the unknown time-dependent surface heat flux in a living skin tissue using the conjugate gradient method. Das et al. (2013) predicted the size and the location of the tumor in a one-dimensional tissue by inverse analysis. Similar studies were also conducted by Das et al. (2008), Partridge and Wrobel (2007), Das and Mishra (2014) and Das and Mishra (2015). Other investigations focused on the simultaneous estimation of various thermo-physical properties of tissues using inverse methods (Luna et al. 2012; Yue et al. 2007; Flach and Özişik, 1989; Huang and Huang, 2004; Jalali et al. 2014). The Pennes model of bioheat transfer equation, which describes the heat transfer in living biological tissues, is widely adopted in variety of biological researches (Trucu et al. 2008, 2009, 2010a, 2010b, 2011; Dehghan and Sabouri, 2012; Fasano et al. 2010; Jasiński et al. 2016). Another mathematical problem that arises during the cryosurgery is the phase change phenomena of the tissue. Indeed, the prediction of the thermo-physical properties during the phase change process by inverse method is also a challenging task (Franquet et al. 2012a, 2012b; Loulou et al. 1999a, 1999b). In all the aforementioned studies however, the prediction of the thermo-physical properties of the biological tissues was not carried out during cryosurgery. The present paper remedies this shortcoming. Its objective is to pursue the previous studies by developing an inverse method capable of estimating the unknown thermo-physical properties during cryosurgery. Once the unknown properties are determined, the isotherms and the moving freezing front are predicted. As a result, the cryosurgical procedure may be monitored and controlled. This paper is organized as follows. First, a Finite-Volume Method (FVM) based on the Pennes bioheat equation and the phase change enthalpy method is developed. The FVM is then validated with an analytical solution available in the open literature. Next, an inverse technique based on the Levenberg-Marquardt Method (LMM) is presented. The (LMM) is combined to the Broyden Method (BM) in order to speed the computations. Finally, the inverse method is thoroughly tested for several freezing scenarios. The effect of the most 5

influential parameters on the inverse predictions and on the course of the cryosurgery is investigated. These parameters are (1) the initial guesses for the unknown LevenbergMarquardt polynomial parameters; (2) the temperature of the cryoprobe; (3) the heat transfer coefficient of the impinging jet of liquid nitrogen; and (4) the noise on the temperature data recorded by the thermocouple probe.

2. Problem statement and assumptions Figure 1 illustrates a schematic of the system under study. The thickness of the diseased tissue layer is

. A cryoprobe is pressed against the layer of tissue (the left boundary at x = 0).

At time t=0, the cryogenic probe temperature is suddenly set equal to the sub-freezing temperature Tc. As a result, heat is transferred from the diseased tissue to the cryoprobe. The temperature of the tissue decreases rapidly and starts freezing on the surface of the cryoprobe. As time passes, the thickness of the frozen layer s(t) increases. The freezing procedure continues until the entire layer of diseased tissue has been deeply frozen and its tumor cells destroyed. The prediction of the heat transfer and of the freezing processes that take place inside the layer of tissue is carried out with a mathematical model that is based on the following assumptions: 

Heat is predominantly transferred from the tissue to the cryoprobe. As a result, a onedimensional heat transfer model is retained (Zhang and Liu, 2002).



The heat transfer inside the tissue is conduction-dominated (Budman et al. 1995; Chua et al. 2007).



The phase change process is non-isothermal. It is characterized by three zones: a frozen solid zone, a mushy zone and an unfrozen zone (Figure 1).



The blood perfusion and the metabolic heat generation are uniform through time and space. These phenomena vanish, however, in the frozen and in the mushy regions.



The blood temperature is set equal to Tb=37 °C.



The heat transfer between the tissue and the blood takes place in the capillaries. The effect of large blood vessels is neglected (Wissler, 1998; Horng et al. 2007).



The liquidus temperature Tliq (the upper limit) and the solidus temperature Tsol (the lower limit) of the biological tissue are set equal to (-1 °C) and (-8 °C), respectively, (Chua and Chou, 2009). 6

3. The direct problem In the direct problem, the thermo-physical properties, the initial and the boundary conditions are all known. The mathematical model for the heat transfer in the tissue layer is given by the following Pennes bio-heat equation (Pennes, 1948):

c

f s T   T   k   b cb wb Tb  T   Qm   H t x  x  t

(1)

Where ρ, c and k are the density, the specific heat and the thermal conductivity of the tissue, respectively. The subscript ‘b’ refers to the blood characteristics. wb and Qm represent the blood perfusion rate and the metabolic heat generation in the tissue, respectively. T is the tissue temperature and Tb is the arterial blood temperature. The last term on the right-hand side of Eq. (1) represents the solid/liquid phase change of the biological tissue during cryosurgery. δH and fs are the enthalpy and the solid fraction during phase change, respectively. The enthalpy δH is defined as  H   (csolid  cliquid )T   L . The solid fraction may be expressed as

0   T  Tliq f s  F T     Tsol  Tliq 1 

 Unfrozen zone 

if

T  Tliq

if

Tsol  T  Tliq  Mushy zone 

if

T  Tsol

(2)

 Frozen solid zone 

At each time-step, the solid fraction fs is updated iteratively as follows:

l



 dF  l 1 1 l f sl 1  f sl    T  F  fs   dT 



(3)

Where F 1 is the inverse function of F. The boundary conditions in Figure 1 are set to

T  x  Ltissue , t   0, x

t

0

(4)

T  x  0, t   Tc ,

t

0

(5) 7

The initial condition is given by

T  x, t=0  37 0C,

x [0, Ltissue ].

(6)

Eq. (1) is solved with a finite-volume method. An implicit time-discretization scheme is adopted. The resulting set of algebraic equations is solved with a Tri-Diagonal-MatrixAlgorithm (TDMA) (Patankar, 1980).

3.1.Validation of the direct model The accuracy of the above mathematical model was first checked with Neumann’s analytical solution for the semi-infinite freezing problem (Carslaw and Jaeger, 1959). The Neumann’s problem was simulated numerically by setting the heat generation and the blood perfusion rate equal to zero in Eq. 1 (Qm=wb=0). The initial temperature was set equal to T0=37 °C. The temperature of the freezing point Tf was set equal to the solidus temperature Tsol (the lower limit). At the left boundary (x=0 m), the temperature of the cryoprobe was imposed Tc. The analytical solution for the time-varying freezing front s(t), provided by Neumann, is:

s(t )  2  f t

(7)

The parameter λ is obtained from the solution of the following transcendental equation

k 1   u 2 k exp     erf     f

where  

 T0  T f   T  T c  f

  exp   2 2   L       erfc    c f T f  Tc 

(8)

f .  f and  u are the thermal diffusivities in the frozen and unfrozen regions, u

respectively. L is the latent heat of fusion. erf and erfc are the error functions. The length of region  0  x  Ltissue  was set equal to Ltissue=0.05 m. The numerical simulation was carried out with a grid made of uniform space increments

=3.125*10-4 m (160 control

8

volumes) and for a constant time-step

. Simulations conducted with finer grid sizes

and smaller time-steps yielded the same results. As an example, Figure 2 depicts the time-varying freezing front for different temperatures of the cryoprobe Tc. It is seen that the agreement between the numerical predictions achieved with the direct model and the analytical solution remains, for all cases, excellent.

4. Inverse model Here, the blood perfusion rate wb, and the thermal conductivities of the frozen and of the unfrozen layers of tissue kf and ku are all unknown. These properties are the most influential as they dictate the heat transfer phenomenon. They are estimated with the inverse model. Once they have been obtained, they are fed to the direct model (section 3 above) that calculates the time-varying temperature distribution T(x,t) and predicts the motion of the freezing front s(t). So the objective of the inverse model is to determine the unknown constant thermal properties of the biological tissue, i.e., wb, kf and ku. To achieve this task, a thermocouple probe is inserted into the diseased tissue at the position x=Ltissue /2=0.025m (Figure 3). The probe records the time-varying temperatures Y  ti  . This information is fed to the inverse model that

 

minimizes the following least squares norm  P

 

I





 P   Y  ti   Tˆ ti , P    i 1



2

(9)



Where Tˆ ti , P is the vector of the estimated temperatures obtained with the inverse model. P = (wb, kf, ku) is the set of the unknown thermal parameters. I is the total number of recorded temperatures. The temperatures recorded by the thermocouple probe, Y  ti  , are called the ‘direct temperatures’. These temperatures are ‘generated’ from the solution of the direct problem for which the parameters P = (wb, kf, ku) are completely specified. Later, these temperatures will be ‘contaminated’ with random noise so as to mimic the effect of measurement errors (see section 5.4).

9

The Levenberg–Marquardt method is adopted for the solution of the present inverse problem (Huang and Huang, 2007; Iljaž and Škerget, 2014). In this iterative method, the incremental value of the unknown thermal parameters ΔP, is expressed as

 

P   J k 

T

J k   k k  

1

 J  Y  T  P   k

T

k

(10)

The superscripts " " and " " refer to the matrix and the vector notation, respectively. The superscript "T" denotes the transpose of the matrix.  is a positive scalar, named ‘Damping k

parameter’.

The

 

k  diag  J k 

T

choice

and

the

update

of

this

parameter

are

discussed

in.

J k  is a diagonal matrix. J k is the Jacobian matrix, also called the 

‘sensitivity matrix’ . It is written as

 T1  w  b  T2  J P   wb    TI  w  b

 

T1 k f T2 k f TI k f

T1  ku   T2   ku    TI  ku 

(11)

The Jacobian matrix plays a crucial role in the estimation of the thermal parameters. There are several approaches for computing the Jacobian (Ozisik and Orlande, 2000). In the present study, the Jacobian matrix is expressed with a finite difference approximation



 

ˆ ˆ Tˆi T ti ; P1 ,..., Pj   Pj  ,...PN  T ti ; P1 ,..., Pj   Pj  ,...PN J ij   Pj 2  Pj 



(12)

10

The subscripts i and j represent the time and the parameter, respectively. The parameter





perturbation  Pj  is set to  1  Pj , where  is a small number. At each iteration, the approximation of the sensitivity coefficients of the Jacobian matrix

J R

I *N



where : I  50 and N  3 requires the solution of the direct problem six times. As

a result, the computation of the Jacobian matrix using a finite difference approximation may become time-consuming. In order to reduce the computational effort, the Jacobian matrix is updated using the Broyden method (Broyden, 1965):

J k  J k 1

Tˆ  Tˆ   J    k

k 1

k 1



Pk 1 PkT1

PkT1Pk 1

(13)

Where J k and J k 1 are the Jacobian matrices at the current and at the previous iteration, respectively, and Pk 1 is the incremental value of the unknown parameters. Further details on this calculation method are provided in (Hafid and Lacroix, 2016, 2017). The algorithms for the direct model and for the inverse model are summarized in Figure 4.

5. Results and discussion In order to monitor the motion of the freezing front during cryosurgery, the blood perfusion rate wb and the thermal conductivities in the frozen and the unfrozen tissues, kf and ku, are first estimated with the inverse model of section 4. This information is then passed to the direct model (section 3 above) that calculates the temperature distribution T(x,t) and predicts the time-varying position of the freezing front s(t). The following numerical simulations were all conducted for a 0.05 m thick biological tissue. The temperature of the cryoprobe was fixed at Tc= -196 °C, and the duration of cryosurgery was limited to 500 s. The physical properties of the biological tissue retained for the simulations are summarized in Table 1 (Deng and Liu, 2005; Kotova et al. 2016) Numerical simulations were carried out for a grid space of

=2.5*10-4 m and a time step of

dt=10 s. Additional simulations were also conducted to ensure that the predictions are grid space and time step independent.

11

All simulations were performed with the Matlab software running on an Intel ® Core(TM) i52520M CPU @ 2,50GHz. The CPU time recorded for a run of the inverse procedure never exceeded 30 s.

To characterize the accuracy of the proposed inverse method, three errors are defined and computed: ErrorP %  100 

ErrorT %  100 

Pexact  Pinverse

(14)

Pexact Texact  Tinverse

(15)

Texact

1 I  s  ti exact  s  ti inverse  RRMSE st  %  100    I i 1  s  ti exact 

2

(16)

P and T stand for the unknown thermal parameter and temperature, respectively. The subscript ‘exact’ refers to the solution obtained from the direct model. The subscript ‘inverse’ refers to the solution estimated with the inverse model. I is the total number of temperature data recorded by the thermocouple probe during cryosurgery. In the present study, this number was set equal to I = 50. 5.1.Effect of the initial guesses This section examines the effect of the initial guesses for the thermal properties Pini   wb ; k f ; ku  on the convergence and the accuracy of the inverse predictions. The

initial guesses were taken as 1/3, 1/5 and 1/7 of that of the exact parameter value (i.e., case#1, case#2 and case#3, respectively). The results of the calculations are shown in Figure 5. They are also summarized in Table 2. It is seen that the inverse estimations are good for all cases. The largest discrepancy occurs in the estimation of the blood perfusion rate wb (4.8%). This is due to its small value with respect to the other thermal properties (kf and ku) considered here (Sawaf and Özisik, 1995). For the remainder of this paper, the initial parameter guesses of case#2 will be retained in all simulations. 12

The inverse predictions for the temperature distribution and for the corresponding ErrorT are illustrated in Figure 6. It is observed that ErrorT remains smaller than 1.5% during the entire duration of the cryosurgical procedure. The maximum value of ErrorT occurs in the mushy zone, that is, in the region where Eq. (1) is nonlinear. The freezing fronts predicted from the direct model and estimated with the inverse model are compared in Figure 7. Examination of this figure reveals that both solutions are in perfect agreement. The corresponding temperature contours calculated with the inverse model are depicted in Figure 8. 5.2. Effect of the freezing temperature The effect of the freezing temperature Tc on the inverse predictions was also checked. Table 3 and Figure 9 show that the inverse method remains stable and accurate in all cases. The temperature contours predicted with the inverse model for Tc= -140 and -230 °C are also illustrated in Figure 10. As expected, the size of the frozen region increases as the freezing temperature Tc drops. 5.3.Effect of heat transfer coefficient h on the inverse estimations For external tumors such as ear and skin tumors, and for cysts, the cryogenic probe is not in direct contact with the surface of the tissue. The surface of the tissue is maintained at a constant sub-freezing temperature Tc by resorting to an impinging jet of liquid nitrogen. In such cases, the temperature field across the layer of tissue and the time-varying location of the freezing front depend on the heat transfer coefficient of the impinging jet (Okajima et al. 2014). The effect of the external heat transfer coefficient on the predictions of the inverse model was examined by replacing Eq. (5) with the following boundary condition (Figure 1):

k

T  x  0, t   h T  0, t   T  x

(17)

T∞ is the temperature of liquid nitrogen. It is fixed at T∞=-196 °C (Ravikumar et al. 1991). Figure 11 illustrates the time evolution of the freezing front, predicted with both the direct and the inverse models, for different values of the heat transfer coefficient h. Once again, the 13

agreement between both models is excellent. The relative root-mean-square error for the freezing front RRMSEs(t) remains, in all cases, below 1%. It is also seen, as expected, that the size of the frozen region increases as the heat transfer coefficient augments. The corresponding temperature contours predicted with the inverse model for h=250 and h=3000 W/m2 °C are depicted in Figure 12. 5.4.Effect of noise on the inverse predictions In order to assess the accuracy of the inverse solution and to obtain the confidence bounds, a statistical analysis for the unknown parameters was performed. The temperatures recorded by the thermocouple probe may be contaminated with measurement errors. For distributed measurement errors with zero mean and a constant variance  2 , the standard deviation of the estimated parameters can be expressed as

 T T  T  Pˆ   diag    T j  P   P

    

1

(18)

Assuming a normal distribution for the errors on the temperatures recorded and 99% confidence, the intervals for the computed parameters Pj are given as

Probability:

 Pˆ  2.576    P j

Pˆj

j ,exact



 Pˆj  2.576  Pˆ

j

  99%

(19)

Pˆj are the estimated values of the unknown polynomial thermal parameters, Pj ,exact , for (j=1,

2, 3), and  Pˆ are the standard deviations given by Eq. (18). j

Let us examine the case with noisy experimental data. The temperatures recorded by the thermocouple probe were contaminated with ‘measurement errors’ (or noise) using the random number generator (randn) i . This measurement error is added to the recorded temperature, which, in the present study, is the temperature provided by the direct model Texact

T (t i )  Texact (t i )  i

(20)

14

Where σ is the standard deviation of the measurement errors. The value for σ is arbitrarily set equal to 1%Tmax and to 2%Tmax.

Tmax is the maximum temperature measured by the

thermocouple probe. Figure 13 compares the measured noisy temperatures (the temperatures generated with the direct model that are contaminated with noise of 2%Tmax) to the estimated temperatures (the temperatures predicted from the inverse model). The confidence bounds of 2.576  Pˆ are also j

shown. Once again, and in spite of the noisy data, the predictions of the inverse model remain accurate and reliable. The effect of noise on the predictions of the freezing front is depicted in Figure 14. As expected, when the noise level rises from 1%Tmax to 2%Tmax, the relative root-mean-square error for the phase change interface RRMSEs(t) slightly increases (from 0.25% to 0.32%). The computational time for noisy data also augments. Nonetheless, the inverse method remains stable and relatively accurate in all cases. 6. Conclusion This article presented an inverse heat transfer method for monitoring the motion of the freezing front from a cryoprobe. With the help of a thermocouple inserted into the layer of diseased tissue, the inverse heat transfer method is able to estimate simultaneously the blood perfusion rate and the thermal conductivities of both frozen and unfrozen tissues. This information is then fed to the Pennes bioheat equation that calculates the time-varying temperature distribution inside the layer of tissue and predicts the motion of the freezing front. The effect of the most influential parameters on the inverse predictions was investigated. These parameters are (i) the initial guesses for the unknown Levenberg-Marquardt polynomial parameters of the thermo-physical properties; (ii) the temperature of the cryoprobe; (iii) the heat transfer coefficient of the impinging jet of liquid nitrogen; and (iv) the noise on the temperature data recorded by the thermocouple probe. Results have shown that the proposed inverse method is a promising alternative to ultrasound and MRI for monitoring the motion of the freezing front during cryosurgery. For all cryogenic scenarios investigated, the predictions of the inverse remained accurate and stable. Future research will focused on the extension of the present 1D-model to multi-dimensional so as to predict the freezing of complex-shaped tumors. Also, in order to monitor in real time the motion of the freezing front, an on-line inverse prediction method based on the Kalman filter will be developed.

15

ACKNOWLEDGEMENTS The authors are grateful to the Natural Sciences and Engineering Research Council of Canada (NSERC) for the financial support.

REFERENCES Broyden C. G., 1965. A class of methods for solving nonlinear simultaneous equations. Mathematics of computation, 19(92), pp. 577-593. Budman H., Shitzer A. & Dayan J., 1995. Analysis of the inverse problem of freezing and thawing of a binary solution during cryosurgical processes. Journal of biomechanical engineering, 117(2), pp. 193-202. Burger M. & Breukelen F. V., 2013. Construction of a low cost and highly sensitive direct heat calorimeter suitable for estimating metabolic rate in small animals. Journal of Thermal Biology, 38(8), pp. 508-512. Cao K. & Lesnic D., 2016. Reconstruction of the perfusion coefficient from temperature measurements using the conjugate gradient method. International Journal of Computer Mathematics, pp. 1-19. Carslaw H. S. & Jaeger J. C., 1959. Conduction of heat in solids. 2 ed. Oxford: Clarendon Press. Carin M. & Jaeger M., 2001. Numerical simulation of the interaction of biological cells with an ice front during freezing. The European Physical Journal of Applied Physics, 16(3), pp. 231-238. Chua K. J. & Chou S. K., 2009. On the study of the freeze–thaw thermal process of a biological system. Applied Thermal Engineering, 29(17), pp. 3696-3709. Cooper T. E. & Trezek G. J., 1971. Correlation of thermal properties of some human tissue with water content. Aerospace medicine, 42(1), pp. 24-27. Chua K. J., Chou S. K. & Ho J. C., 2007. An analytical study on the thermal effects of cryosurgery on selective cell destruction. Journal of biomechanics, 40(1), pp. 100-116. Das K., Singh R. & Mishra S. C., 2013. Numerical analysis for determination of the presence of a tumor and estimation of its size and location in a tissue. Journal of thermal biology, 38(1), pp. 32-40. Das K. & Mishra S. C., 2014. Non-invasive estimation of size and location of a tumor in a human breast using a curve fitting technique. International Communications in Heat and Mass Transfer, Volume 56, pp. 63-70.

16

Das K. & Mishra S. C., 2015. Simultaneous estimation of size, radial and angular locations of a malignant tumor in a 3-D human breast–A numerical study. Journal of thermal biology, Volume 52, pp. 147-156. Das R., Mishra S. C. & Uppaluri R., 2008. Multiparameter estimation in a transient conduction-radiation problem using the lattice Boltzmann method and the finite-volume method coupled with the genetic algorithms. Numerical Heat Transfer, 53 (12), pp. 13211338. Dehghan M. & Sabouri M., 2012. A spectral element method for solving the Pennes bioheat transfer equation by using triangular and quadrilateral elements. Applied mathematical modelling, 36(12), pp. 6031-6049. Delhomme G., et al., 1992. Thermal diffusion probes for tissue blood flow measurements. Sensors and Actuators B: Chemical, 6 (1-3), pp. 87-90. Deng Z. S. & Liu J., 2005. Numerical simulation of selective freezing of target biological tissues following injection of solutions with specific thermal properties. Cryobiology, 50(2), pp. 183-192. Fasano A., Hömberg D. & Naumov D., 2010. On a mathematical model for laser-induced thermotherapy. Applied Mathematical Modelling, 34(12), pp. 3831-3840. Flach G. P. & Özişik M. N., 1989. Inverse heat conduction problem of simultaneously estimating spatially varying thermal conductivity and heat capacity per unit volume. Numerical Heat Transfer, 16(2), pp. 249-266. Franquet E., Bedecarrats J.P., Haillot D. & Dumas J. P., 2012a. Determination of the enthalpy of phase change materials by inverse method from calorimetric experiments. Applications to pure substances or binary solutions. Journal of Physics: Conference Series, p. 012135. Franquet E., Gibout S., Bedecarrats J.P., Haillot D. & Dumas J. P., 2012b. Inverse method for the identification of the enthalpy of phase change materials from calorimetry experiments. Thermochimica acta, Volume 546, pp. 61-80. Gage A. A. & Baust J. G., 2007. Cryosurgery for tumors. Journal of the American College of Surgeons, 205(2), pp. 342-356. Gosselin L., Tye-Gingras M. & Mathieu-Potvin F., 2009. Review of utilization of genetic algorithms in heat transfer problems. International Journal of Heat and Mass Transfer, 52(9), pp. 2169-2188. Grabski J. K., Lesnic D. & Johansson B. T., 2016. Identification of a time-dependent bio-heat blood perfusion coefficient. International Communications in Heat and Mass Transfer, Volume 75, pp. 218-222. Hafid M. & Lacroix M., 2016. An inverse heat transfer method for predicting the thermal characteristics of a molten material reactor. Applied Thermal Engineering, Volume 108, p. 140–149. Hafid M. & Lacroix M., 2017. Inverse heat transfer prediction of the state of the brick wall of a melting furnace. Applied Thermal Engineering, Volume 110, pp. 265-274. 17

Hazanee A. & Lesnic D., 2014. Determination of a time-dependent coefficient in the bioheat equation. International Journal of Mechanical Sciences, Volume 88 , pp. 259-266. Huang C. H. & Huang C. Y., 2004. An inverse biotechnology problem in estimating the optical diffusion and absorption coefficients of tissue. International journal of heat and mass transfer, 47(3), pp. 447-457. Huang C. H. & Huang C. Y., 2007. An inverse problem in estimating simultaneously the effective thermal conductivity and volumetric heat capacity of biological tissue. Applied mathematical modelling, 31(9), pp. 1785-1797. Huber P. E., et al., 2001. A new noninvasive approach in breast cancer therapy using magnetic resonance imaging-guided focused ultrasound surgery. Cancer research, 61(23), pp. 8441-8447. Huntul M. J., Lesnic D. & Hussein M. S., 2017. Reconstruction of time-dependent coefficients from heat moments. Applied Mathematics and Computation, Volume 301, pp. 233-253. Horng T. L., et al., 2007. Effects of pulsatile blood flow in large vessels on thermal dose distribution during thermal therapy. Medical physics, 34(4), pp. 1312-1320. Iljaž J. & Škerget L., 2014. Blood perfusion estimation in heterogeneous tissue using BEM based algorithm. Engineering Analysis with Boundary Elements, Volume 39, pp. 75-87. Jaeger M. & Carin M., 2002. The front-tracking ALE method: application to a model of the freezing of cell suspensions. Journal of Computational Physics, 179(2), pp. 704-735. Jalali A., Ayani M. B. & Baghban M., 2014. Simultaneous estimation of controllable parameters in a living tissue during thermal therapy. Journal of thermal biology, Volume 45, pp. 37-42. Jasiński M., Majchrzak E. & Turchan L., 2016. Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag equation. Applied Mathematical Modelling, 40(2), pp. 750-762. Khanday M. A. & Hussain F., 2015. Explicit formula of finite difference method to estimate human peripheral tissue temperatures during exposure to severe cold stress. Journal of thermal biology, Volume 48, pp. 51-55. Korpan N. N., 2012. Basics of cryosurgery. s.l.:Springer Science & Business Media. Kotova Т. G., et al., 2016. Calculation of Effective Freezing Time in Lung Cancer Cryosurgery Based on Godunov Simulation. Современные технологии в медицине, 8(1), pp. 48-53. Lee C. Y. & Bastacky J., 1995. Comparative mathematical analyses of freezing in lung and solid tissue. Cryobiology, 32(4), pp. 299-305. Lesnic D., 2009. Identification of the time-dependent perfusion coefficient in the bio-heat conduction equation. Journal of Inverse and Ill-posed Problems, 17(8), pp. 753-764. Lesnic D. & Ivanchov M., 2015. Determination of the time-dependent perfusion coefficient in the bio-heat equation. Applied Mathematics Letters, Volume 39, pp. 96-100. 18

Lee H. L., Lai T. H., Chen W. L., & Yang Y. C., 2013. An inverse hyperbolic heat conduction problem in estimating surface heat flux of a living skin tissue. Applied Mathematical Modelling, 37(5), pp. 2630-2643. Luna J. M., Romero-Mendez R., Hernandez-Guerrero A. & Elizalde-Blancas F., 2012. Procedure to estimate thermophysical and geometrical parameters of embedded cancerous lesions using thermography. Journal of biomechanical engineering, 134(3), p. 031008. Loulou T., Artyukhin E. A. & Bardon J. P., 1999a. Estimation of thermal contact resistance during the first stages of metal solidification process: I—experiment principle and modelisation. International Journal of Heat and Mass Transfer, 42(12), pp. 2119-2127. Loulou T., Artyukhin E. A. & Bardon J. P., 1999b. Estimation of thermal contract resistance during the first stages of metal solidification process: II—experimental setup and results. International journal of heat and mass transfer, 42(12), pp. 2129-2142. Marquardt D. W., 1963. An algorithm for least-squares estimation of nonlinear parameters. Journal of the society for Industrial and Applied Mathematics, 11 (2), pp. 431-441. Malawer M. M., et al., 1999. Cryosurgery in the Treatment of Giant Cell Tumor: A Long term Followup Study. Clinical orthopaedics and related research, Volume 359, pp. 176-188. Okajima J., Komiya A. & Maruyama S., 2014. 24-gauge ultrafine cryoprobe with diameter of 550μm and its cooling performance. Cryobiology, 69(3), pp. 411-418. Ozisik M. N. & Orlande H. R. B., 2000. Inverse Heat Transfer. New York: Taylor and Francis. Partridge P. W. & Wrobel L. C., 2007. An inverse geometry problem for the localisation of skin tumours by thermal analysis. Engineering Analysis with Boundary Elements, 31(10), pp. 803-811. Patankar S., 1980. Numerical heat transfer and fluid flow. Florida: CRC press. Pennes H. H., 1948. Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of applied physiology, 1(2), pp. 93-122. Pfleiderer S. O., et al., 2002. Cryotherapy of breast cancer under ultrasound guidance: initial results and limitations. European radiology, 12(12), pp. 3009-3014. Ravikumar T. S., Steele G., Kane R. & King V., 1991. Experimental and clinical observations on hepatic cryosurgery for colorectal metastases. Cancer research, 51(23), pp. 6323-6327. Sawaf B. & Özisik M. N., 1995. Determining the constant thermal conductivities of orthotropic materials by inverse analysis. Int. communications in heat and mass transfer, 22(2), pp. 201-211. Trucu D., Ingham D. B. & Lesnic D., 2008. Inverse space-dependent perfusion coefficient identification. Journal of Physics: Conference Series. IOP Publishing, 135(1), p. 012098. Trucu D., Ingham D. B. & Lesnic D., 2009. An inverse coefficient identification problem for the bio-heat equation. Inverse Problems in Science and Engineering, 17(1), pp. 65-83.

19

Trucu D., Ingham D. B. & Lesnic D., 2010a. Inverse temperature-dependent perfusion coefficient reconstruction. International Journal of Non-Linear Mechanics, 45(5), pp. 542549. Trucu D., Ingham D. B. & Lesnic D., 2010b. Space-dependent perfusion coefficient identification in the transient bio-heat equation. Journal of Engineering Mathematics, 67(4), pp. 307-315. Trucu D., Ingham D. B. & Lesnic D., 2011. Reconstruction of the space-and time-dependent blood perfusion coefficient in bio-heat transfer. Heat Transfer Engineering, 32(9), pp. 800810. Wissler E. H., 1998. Pennes’ 1948 paper revisited. Journal of Applied Physiology, 85 (1), pp. 35-41. Yue K., Zhang X. & Yu F., 2007. Simultaneous estimation of thermal properties of living tissue using noninvasive method. International Journal of Thermophysics, 28(5), pp. 14701489. Zhang Y. T. & Liu J., 2002. Numerical study on three-region thawing problem during cryosurgical re-warming. Medical engineering & physics, 24(4), pp. 265-277. Zhao X. & Chua K. J., 2012. Studying the thermal effects of a clinically-extracted vascular tissue during cryo-freezing. Journal of Thermal Biology, 37(8), pp. 556-563.

20

Table 1. Properties of the biological tissue.

Thermal conductivity of unfrozen tissue Thermal conductivity of frozen tissue Specific heat of unfrozen tissue Specific heat of frozen tissue Density of unfrozen tissue Density of frozen tissue Latent heat of freezing Lower phase transition temperature Upper phase transition temperature Blood perfusion rate Metabolic heat generation Specific heat of blood Density of blood Blood temperature

Symbol ku kf cu cf ρu ρf L Tsol Tliq wb Qm cb ρb Tb

Value 0.5 2 3600 1800 1000 1000 250 −8 −1 0.0005 4200 3640 1000 37

Unit W/m °C W/m °C J/kg °C J/kg °C kg/m3 kg/m3 kJ/kg °C °C ml/s ml W/m3 J/kg °C kg/m3 °C

21

Table 2. Effect of the initial guesses on the predictions of the thermal parameters. Case#1 Pinitial= Pexact/3

Case#2 Pinitial= Pexact/5

Case#3 Pinitial= Pexact/7

Pexact

Pinverse

ErrorP%

Pinverse

ErrorP%

Pinverse

ErrorP

wb (ml/ml*s)

5*10-4

4.92*10-4

1.6

4.96*10-4

0.8

4.76*10-4

4.8

kf (W/m °C)

2

1.9998

0.01

1.9999

0.005

1.9992

0.04

ku (W/m °C)

0.5

0.5

0

0.4999

0.02

0.4997

0.06

22

Table 3. Effect of the freezing temperature on the parameter estimation.

Tc=-140 °C

Tc=-196 °C

Tc=-230 °C

PExact

PInverse

Error%

PInverse

Error%

PInverse

Error%

wb

5*10-4

4.93*10-4

1.4

4.96*10-4

0.8

4.97*10-4

0.6

kf

2.00

1.9998

0.01

1.9999

0.005

2

0

ku

0.50

0.4999

0.02

0.4999

0.02

0.4999

0.02

23

Figure 1. Schematic of the system. The time-varying position of the freezing front is s(t).

24

0.03 Numerical solution Analytical solution

Interface position (m)

0.025

0.02

0.015 o

Tc=-120 ( C) o

Tc=-160 ( C)

0.01

o

Tc=-196 ( C) o

Tc=-240 ( C) 0.005

0

0

50

100

150

200

250

300

350

400

450

500

Time (s)

Figure 2. Time-varying freezing front position for different temperatures of the cryoprobe.

25

Cryoprobe

Tsol

Tliq

Frozen tissue

kf=?

Tc

Thermocouple probe (x=Ltissue /2)

.

Mushy region

Unfrozen tissue

ku=? wb=?

Ltissue x

s(t)=?

Figure 3. The inverse problem. wb , ku and kf are unknown. They are estimated from the temperatures recorded by thermocouple probe.

26



Thermo-physical properties Pexact  wb ; k f ; k u



Setup initial temperature T and solid fraction fs Store old values: T and fs t=t+dt Solve bio-heat equation

Update fs

t
T(x,t), s(t)

Initial guesses Pinitial

Thermo-physical properties Pexact

P k 1

Noisy data

Direct Model

Sensor

Y  ti 



Tˆ ti , P

Pfit  wb ; k f ; ku 

Interface position s(t) and temperature distribution T(x,t)

Objective function   P 



Minimization?

yes

Convergence criteria ? no ∂Ti/∂Pj Full computational of J ?

yes

no Broyden update of J

Figure 4. Direct and inverse heat transfer algorithms. 27

2

x 10

-3

Case#1 Case#2 Case#3

1.5

wb (ml/ml. s)

1 0.5 0 -0.5 -1 -1.5 -2

0

2

4

6

8

10

12

14

16

18

14

16

18

Iteration number

(5-a) 2.2 2 1.8

k f (W /m °C)

1.6

Case#1 Case#2 Case#3

1.4 1.2 1 0.8 0.6 0.4 0.2

0

2

4

6

8

10

12

Iteration number

(5-b)

28

1.8 Case#1 Case#2 Case#3

1.6 1.4

k u (W /m °C)

1.2 1 0.8 0.6 0.4 0.2 0

0

2

4

6

8

10

12

14

16

18

Iteration number

(5-c) Figure 5. Convergence of the unknown thermal parameters.

29

Inverse temperature (K)

50 0 0 -50 -50 -100 -100 -150 -150 500 400

0.05 300

Time (s)

0.04 0.03

200

0.02

100

0.01

x (m)

0

(a)

ErrorT (%)

1.5

1

0.5

0 500 400

0.05 300

Time (s)

0.04 0.03

200

0.02

100

0.01 0

x (m)

(b)

Figure 6. Inverse predictions of the temperature distribution (a) and ErrorT (b). 30

0.03 Exact solution inverse solution

Interface position (m)

0.025

0.02

0.015

0.01

0.005

0

0

50

100

150

200

250

300

350

400

450

500

Time (s)

Figure 7. Predicted freezing front without noise.

31

40

40

0.045

20 0

0.04

-20

40

0.035

-40

20 0

0.025

20 0

0.02 20 0

0.015 0.01 200 ---468000 0 ---1 1 246000 -1 -180 100

0 20

x (m)

0.03

0.005

-2 0 -4 0 -6 0 -8 0 -1 00 -1 20 -1 40 -160 -180

200

300

-60 -80 -100

-2 0 -4 0 -6 0 -8 0 -1 00 -1 20 -140 -160 -180

400

-120 -140 -160 -180

500

Time (s) Figure 8. Temperature contours predicted with the inverse model.

32

0.03 Exact solution inverse solution

Interface position (m)

0.025

0.02 Tc=-230 °C

0.015 Tc=-196 °C

0.01

Tc=-140 °C

0.005

0

0

50

100

150

200

250

300

350

400

450

500

Time (s) Figure 9. Motion of the freezing front for different temperatures Tc.

33

40

0.045

20

0

0.04

0.04

0

0.035

40

0.035

40

-20

-50 20 0

0

20 0

0.02 20 0

0.015 0.01 0 -24 00 -6 --8000 0 -1 -12

100

-2 0 -4 0 -6 0

-2 0 -4 0 -6 0 -8 0 -100 -120

200

300

Tc=-140 °C

-60 -80

-80 -100 -120

Time (s)

-40

400

500

0.025

20 0 -2 0 -4 0 -6 0 -8 0 -1 00 -1 20 -140 -160 -180 -200 -220

0.02 20 0 -2 0-4 0 -6 0 -8 0 0 -1 0 -1 20 -1 40 -1 60 -180 -200 -220

0.015

-100

0.01

-120

0.005

0 20

20

0.025

x (m)

0.03

0 20

x (m)

0.03

0.005

40

40

0.045

0 -4 00 --81620000 4600 1 ---1 -1-20 8 0

100

200

300

400

-2 0

-100

-150

-200

500

Time (s)

Tc=-230 °C

Figure 10. Isotherms predicted with the inverse model.

34

0.03 Exact solution inverse solution

Interface position (m)

0.025

0.02

0.015

2

h=3000 W/m °C 2

0.01

h=750 W/m °C

2

0.005

0

h=250 W/m °C

0

50

100

150

200

250

300

350

400

450

500

Time (s)

Figure 11. Effect of the heat transfer coefficient.

35

h=250 W/m2 °C

h=3000 W/m2 °C

Figure 12. Temperature contours predicted with the inverse model.

36

40 35 20 30

Temperature °C

0

25 20

-20

40

60

80

-40

-60 Temperature measurements (Direct) Temperature estimated (Inverse) Confidence bounds

-80

-100

50

100

150

200

250

300

350

400

450

500

Time [s]

Figure 13. Measured and predicted inverse temperatures.

37

0.03 Exact Inverse  =1% Tmax

Interface position (m)

0.025

Inverse  =2% Tmax

0.02

0.015 0.0163 0.01 0.0162 0.0161

0.005

183 0

0

50

100

150

200

250

184

300

185

350

186

400

450

500

Time (s)

Figure 14. Predicted position of the freezing front with noise.

38

Highlights:



Bioheat transfer in living tissue are considered.



An inverse analysis is proposed to improve the cryosurgery.



The Levenberg-Marquardt-Method combined with the Broyden-Method is adopted.



The prediction of the thermo-physical proprieties during cryosurgery is performed.



The effect of the most influential parameters on the inverse predictions is examined.

39