- Email: [email protected]
Source term prediction in a multilayer tissue during hyperthermia
Author’s Accepted Manuscript Source term prediction in a multilayer tissue during hyperthermia M. Baghban, M.B. Ayani
www.elsevier.com/locate/jtherbio
PII: DOI: Reference:
S0306-4565(15)30065-6 http://dx.doi.org/10.1016/j.jtherbio.2015.07.006 TB1659
To appear in: Journal of Thermal Biology Received date: 2 March 2015 Revised date: 14 July 2015 Accepted date: 14 July 2015 Cite this article as: M. Baghban and M.B. Ayani, Source term prediction in a multilayer tissue during hyperthermia, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2015.07.006 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.
Source term prediction in a multilayer tissue during hyperthermia
M. Baghban, M. B. Ayani Department of Mechanical Engineering Faculty of Engineering Ferdowsi University of Mashhad Mashhad, Iran P.O. Box 91775-1111 Phone No.: 0098 9155117671 Fax: 0098 513 8763304 Email: [email protected]
Abstract: One of the major challenges in the use of hyperthermia to treat cancer is determining the desired heating power of external source in such a way that the thermal injury is confined to the unhealthy tissue. In this study, an inverse method based on the sequential method is proposed to estimate the desired heating power as a function of time for a successful hyperthermia treatment. In order to simulate the measured temperature, the direct problem is solved for a multilayer skin tissue to obtain the temperature data at the skin surface. These data are employed in the inverse problem to estimate the heating power of external source. Two examples are considered to examine the accuracy of the inverse analysis. In addition, the effect of measurement errors is investigated. Results show that the proposed inverse algorithm is able to determine the desired heating power of external source accurately, even in the presence of measurement errors. However, for noisy data, more temperature measurements are required to achieve reliable results.
Keywords: Inverse problem, Sequential method, Hyperthermia, Multilayer skin tissue, Heating power, Heat source.
1. Introduction Hyperthermia is a method of cancer treatment in which the temperature of tumors increased in the range of 42°C–46 °C for a specified period of time to destroy them while the healthy parts are not destroyed (Field and Bleehen, 1979; Overgaard, 1989). In the clinical application, hyperthermia can be used in two different ways, local hyperthermia and regional hyperthermia. In the local hyperthermia, very high temperatures are applied to destroy the tumor while in the regional hyperthermia the temperatures increase a few degrees to help other cancer treatments such as chemotherapy and radiotherapy (Hildebrandt et al., 2002).
Nomenclature
∆x Y Z
specific heat, J/kg K objective function heat transfer coefficient, W/m2 oC thermal conductivity, W/m oC length, m number of nodes at spatial coordinate number of temporal measurements number of future time steps source term, W/m3 time-dependent part of heat source, heating power, W/m2 temperature, K, oC temporal coordinate, s increment of time domain, s spatial coordinate, m increment of spatial coordinate, m measured temperature, K, oC sensitivity matrix
Greek symbols β
scattering coefficient, 1/m
c F(q) h k L M N r S(x,t) q(t) T(x,t) t ∆t
x
ρ ζ σ ω
density random variable standard deviation of the measurement errors blood perfusion
Subscript b es d E f if s P W ∞
blood estimated desired east node final time interface s-th layer central node west node cooling fluid
Superscript n *
time level Initial guess
To have a successful thermal therapy, it is necessary to control the temperature. For this goal, many researchers attempted to solve inverse bio-heat conduction in estimating the boundary condition and heat source. Loulou and Scott (2002) used the conjugate gradient method to optimize the heat source that results in the desired thermal dose. Kuznetsov (2006) applied the minimum principal of Pontryagin and determined the optimal volumetric heating rate. Yang (2011) estimated the boundary condition in hyperbolic bio-heat conduction by using the sequential method. Dhar et al. (2012) investigated analytically optimal distribution of time dependent heating power in a multilayer tissue. Lee et al. (2013) developed an inverse hyperbolic analysis for estimating the unknown time-dependent surface heat flux of a living skin tissue. Yang (2014a) predicted the boundary condition in two-dimensional non-Fourier multilayer bioheat conduction. Lee et al (2014) used the conjugate gradient method to estimate the unknown time-dependent surface heat flux in a multilayer tissue based on the hyperbolic model of bio-heat
conduction. Jalali et al. (2014) used the conjugate gradient method to estimate simultaneously the heat source and heat transfer coefficient at the skin surface for the desired temperature distribution. Yang (2014b) determined the heat strength that destroys tumor tissue but does not damage normal tissue. This study deals with determining the unknown time-dependent heating power of external source using the sequential algorithm of the inverse bio-heat transfer. In order to simulate the measured temperature, the direct problem is solved for a one-dimensional multilayer tissue to obtain the temperature data at the skin surface. The accuracy of the inverse method is checked by utilizing two examples.
2. Direct problem The physical model applied in this work is sketched in Fig. 1. Assume that the thermal properties are constant and the metabolic heat generation is negligible compared to the heat source due to external heating; the thermal model of bio-heat transfer in the s-th layer of the tissue is given as (Pennes, 1948): ρscs
∂T ∂t
= ks
∂2T ∂x 2
+ ωbscbs (Tb −T ) + S (x ,t )
(1)
where ρ, c and k are the density, specific heat and thermal conductivity of the tissues, respectively. ωb, cb and Tb denote the blood perfusion rate, specific heat and the arterial blood temperature, respectively. S(x,t) is the heat source term where heat can be generated from external sources such as laser or ultrasound. The initial condition and boundary conditions are expressed as:
T (x ,t ) = 37 °C for t = 0
(2)
∂T (x , t ) = h (T ∞ −T (x , t ) ) ∂x
−k
(3)
at x = 0
T (x , t ) = 37 °C at x = L
(4)
where T∞ and h represent the temperature and heat-transfer coefficient of the cooling fluid on the skin surface, respectively. At the interface of two layers, the temperature and heat flux are continuous. The boundary conditions at the interface of i-th and j-th layers (x=xif) are as follows:
T i (x , t ) = T j (x , t ) −k i
(5)
at x = x if
∂T j ( x , t ) ∂T i (x , t ) = −k j ∂x ∂x
(6)
at x = x if
The fully implicit finite difference method (Ferzige and Peric, 2002) is used to solve the direct problem, subject to the initial condition and the boundary conditions. Therefore, the governing equation (Eq. (1)) can be expressed as the following: k s − (∆x )2
ρ c 2k s T E n +1 + s s + + ωbs cbs ∆t (∆x )2
ρs c s n T P ∆t
k s T P n +1 − (∆x )2
TW n +1 =
(7)
+ (ωbs cbs )Tb + S (x ,t )
where TPn+1 is the temperature of node P at the time n+1 (see Fig. 2). ∆x and ∆t are the increments of spatial coordinate and time domain, respectively. TE and TW are the temperature of the east and west nodes, respectively. In the direct problem, the source term is known, while in the inverse problem, it is considered as an unknown value and needs to be predicted by measuring the temperatures at the skin surface.
3. Inverse problem
In this study, the sequential method combined with the concept of the future time (Beck et al., 1985) is used to estimate the heat source term. According to the well-known Beer's Law, the heat source term is given by (Deng and Liu, 2002),
S (x ,t ) = β q (t )e −βx
(8)
where β is the scattering coefficient. q(t) is the time-dependent part of source term, which is the heating power of the external source. The objective of the inverse analysis is to estimate q(t) from the knowledge of temperature measurement taken at the skin surface (x=0). To estimate the transient heating power, the heating power is broken up into ‘‘n’’ elements. Therefore, the governing equation at t=tm is rewritten by the following equations: ρscs
∂T (x ,t m ) ∂t
= ks
∂2T (x ,t m ) ∂ 2x
+ ωbscbs (Tb −T (x ,t m )) + qm +1β e −βx
(10)
T (x , t ) = T m −1 for t = t m −1
∂T (x , t m ) = h (T ∞ −T (x , t m ) ) ∂x
−k
T (x , t m ) = 37 °C
(9)
at x = 0
(11)
(12)
at x = L
At the interface of i-th and j-th layers (x=xif), the boundary conditions are as follows:
T i (x , t m ) = T j (x , t m ) −k i
(13)
at x = x if
∂T j (x , t m ) ∂T i ( x , t m ) = −k j ∂x ∂x
at x = x if
(14)
The temperature distribution at t=tm is determined by solving Eqs. (9)-(14). In order to estimate qm+1, it is supposed that the values q1, q2, ..., qm have been previously calculated. In addition, the heating power components at the future time are assumed to be constant over "r" future time steps, i. e,
(15)
qm+2 = ... = qm+r −1 = qm +r = qm +1
The objective function for the identification of qm+1 is expressed by the sum of square residuals between the measured and estimated temperatures as follows: r
F (qm+1) = ∑[Y m+i −Tm+i (qm+1)]
2
(16)
i =1
Equation (16) is minimized with respect to the qm+1 to obtain, r dF = −2∑ [Y m+i −Tm+i (qm+1)][ Z i ] dqm+1 i =1
(17)
where the step sensitivity coefficients, Z, at the sensor location defined as:
∂T ∂Tm+2 ∂Tm+r Z = m+1 … ∂qm+1 ∂qm+1 ∂qm+1
(18)
Employing a Taylor series expansion of Tm+1(qm+1) about an initial guess, q*, qm+1 is calculated as: (19)
r
∑(Y qm +1 = q* + i =1
m +i
−Tm +i )Z i
r
∑Z
2 i
i =1
In linear problems such as the present study, Eq. (19) is an accurate equation for estimating the unknown heating power.
3.1 Sensitivity Problem To calculate the sensitivity coefficients (Eq. (18)), the direct problem given by Eqs. (9)-(14) is differentiated with respect to qm+1 as: ρscs
∂Z (x ,t m ) ∂t
Z (x , t m ) = 0
= ks
∂2Z (x ,t m )
for
∂2x
−ωbscbs Z (x ,t m ) + β e −βx
t m = t m −1
(20)
(21)
k
∂Z (x , t m ) = h Z (x , t m ) at x = 0 ∂x
(22)
Z (x , t m ) = 0 at x = L
(23)
Z i (x , t m ) = Z j (x , t m ) at x = x if
(24)
−k i
∂Z j ( x , t m ) ∂Z i (x , t m ) = −k j ∂x ∂x
at x = x if
(25)
The finite difference method is used to solve the sensitivity problem.
4. Results and discussion The sequential algorithm developed in the present study is tested by using it for the problem of estimating the time-dependent heating power, q(t), in a one-dimensional multilayer tissue. Since experimental data are not available, a real heat source is assumed to calculate the temperatures at the measurement location (x=0). The results are considered as the computed temperature, Y. In order to simulate the measured temperature containing measurement errors, random errors of standard deviation, σ, are added to the temperature computed from the solution of the direct problem. It can be shown in the following equation,
Y =Y ex +σζ
(26)
where ζ is a random variable within -2.576 to 2.576 for a 99% confidence bounds. The physical and thermophysical properties of all layers are listed in Table 1 (Liu and Cheng, 2008). The skin surface is a convection boundary with the convection coefficient h=200 W/m2 oC and the cooling temperature T∞=5 oC. The scattering coefficient for all layers is taken as β=100 m-1 (Dhar et. al, 2012). The final time is chosen as tf =1000s. The direct solution for a one-layer tissue is validated by Deng's analytical solution (Deng and Liu, 2002) by assuming the same conditions. To demonstrate the accuracy and efficiency of the inverse solution in predicting the
time dependent heating power, two different functions involving abrupt changes are considered in this study. In all examples, the estimated time-dependent heating power, qes(t), obtained with an initial guess q*(t)=0.
Example 1. The desiered transient heating power, qd (t), is assumed in the following form,
500 qd (t ) = 1000 500
tf 3 tf 2t ≤t ≤ f 3 3 2t f < t ≤ tf 3 0 ≤t <
(27)
Fig. 3 displays the effect of future time steps number on the accuracy of estimation without measurement errors. It can be seen clearly that increasing the number of future time steps, improves the estimating stability but decreases the accuracy. To determine the accuracy of the inverse solution, the relative percentage error is defined by: N
qes −qd ×100, qd m=1
ERR = ∑
(28)
where N is the number of temporal measurements. The relative percentage errors for different number of future time steps and measurement errors are listed in Table 2. It is obvious from the results that r=24 could present a stable and accurate estimation when σ=0. In addition, it can be concluded that the heat source term can be estimated accurately, even with noisy data. However, for noisy data, more number of future time steps is required to get an accurate solution. The
influence of the measurement errors on the solution for example 1, when the future time steps number is 24, is shown in Fig. 4. For r=24, the relative errors for q(t) are computed as 2.425%, 3.524% and 4.808%, for σ=0, σ=0.1 and σ=0.2,
respectively. The results display that the
proposed solution is an accurate and stable method to determine the time-dependent heating power.
Example 2. The desired time-dependent heating power, qd (t), is taken in the following form,
t 0 ≤t ≤ f 500 4 2000(t − t f )/ t + 500 t f ≤ t ≤ t f f 4 4 2 qd (t ) = −2000(t − 3t f )/ t + 500 t f ≤ t ≤ 3tf f 4 2 4 3t f 500 ≤t ≤ tf 4
(29)
The effects of future time steps and measurement errors on the accuracy and stability of estimation is demonstrated in Figs. 5 and 6, respectively. In addition, the relative percentage errors for different number of future time steps and measurement errors are listed in Table 3. It can be seen, the relative percentage error is decreases to the minimum value and then increases when the future time steps number increases. However, increasing r results in a more stable solution. Also, it is clear that the relative error is increased with the increase of the measurement errors. For r=24, the relative error decreases 49% as the measured error reduces from σ=0.2 to
σ=0. The results show that the estimated values of the time-dependent heating power agree well with the desired values even for the case with the measurement errors.
5. Conclusion
An inverse algorithm based on the sequential method was used to estimate the time-dependent part of source term in a one-dimensional multilayer skin tissue from the temperature measurement at the skin surface. The Pennes bio-heat equation was applied to simulate thermal behavior in the tissue. Two test cases were presented to check the ability and accuracy of the proposed inverse solution. The effect measurement noise on the accuracy of the estimation was studied. The results showed that the agreement between the estimated and desired values of the heat source was excellent, even with noisy data.
References Beck, J. V., Blackwell, B., Clair Jr, C. R. S., 1985. Inverse heat conduction: Ill-posed problems. Wiley, New York. Deng, Z. S., Liu, J., 2002. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. Journal of biomechanical engineering, 124 (6), 638-649. Dhar, P., Dhar, R., Dhar, R., 2012. Analytical study on optimal distribution of heating power in hyperthermia. International Communications in Heat and Mass Transfer, 39 (3), 419-423. Ferziger, J. H., Peric¸ M., 2002. Computational Methods for Fluid Dynamics. Springer, New York. Field, S. B., Bleehen, N. M., 1979. Hyperthermia in the treatment of cancer. Cancer treatment reviews, 6 (2), 63-94. Wust, P., Hildebrandt, B., Sreenivasa, G., Rau, B., Gellermann, J., Riess, H., Felix, R., Schlag, P. M., 2002. Hyperthermia in combined treatment of cancer. The lancet oncology, 3 (8), 487-497.
Jalali, A., Ayani, M. B., Baghban, M., 2014. Simultaneous estimation of controllable parameters in a living tissue during thermal therapy. Journal of thermal biology, 45, 37-42. Kuznetsov, A. V., 2006. Optimization problems for bioheat equation. International Communications in Heat and Mass Transfer, 33 (5), 537-543. 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), 2630-2643.
Lee, H. L., Chen, W. L., Chang, W. J., Yang, Y. C., 2015. Estimation of surface heat flux and temperature distributions in a multilayer tissue based on the hyperbolic model of heat conduction. Computer methods in biomechanics and biomedical engineering, 18 (14), 15251534. Liu, K. C., Cheng, P. J., 2008. Finite propagation of heat transfer in a multilayer tissue. Journal of Thermophysics and Heat Transfer, 22 (4), 775-782. Loulou, T., Scott, E. P., 2002. Thermal dose optimization in hyperthermia treatments by using the conjugate gradient method. Numerical Heat Transfer: Part A: Applications, 42 (7), 661-683. Overgaard, J., 1989. The current and potential role of hyperthermia in radiotherapy. nternational Journal of Radiation Oncology • Biology • Physics, 16 (3), 535-549. Pennes, H. H., 1948. Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of applied physiology, 1 (2), 93-122. Yang, C. Y., 2011. Boundary estimation of hyperbolic bio-heat conduction. International Journal of Heat and Mass Transfer, 54 (11), 2506-2513.
Yang, C. Y., 2014a. Boundary prediction of bio-heat conduction in a two-dimensional multilayer tissue. International Journal of Heat and Mass Transfer, 78, 232-239. Yang,
C.
Y.,
2014b.
Determining
the
heat
strength
required
in
hyperthermia
treatments. International Communications in Heat and Mass Transfer, 57, 282-285.
Figure Captions: Fig. 1. Schematic of a one-dimensional multilayer tissue. Fig. 2. Schematic of a node. Fig. 3. Comparison between the desired and estimated heating power for different future time steps in example 1 without measurement errors (σ=0). Fig. 4. Comparison between the desired and estimated heating power for different measurement errors in example 1 when r=24. Fig. 5. Comparison between the desired and estimated heating power for different future time steps in example 2 without measurement errors (σ=0). Fig. 6. Comparison between the desired and estimated heating power for different measurement errors in example 2 when r=24.
Table Captions: Table 1. Physical and thermophysical properties (Liu and Cheng, 2008). Table 2. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 1. Table 3. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 2.
Fig. 1. Schematic of a one-dimensional multilayer tissue.
Fig. 2. Schematic of a node.
Fig. 3. Comparison between the desired and estimated heating power for different future time steps in example 1 without measurement errors (σ=0).
Fig. 4. Comparison between the desired and estimated heating power for different measurement errors in example 1 when r=24.
Fig. 5. Comparison between the desired and estimated heating power for different future time steps in example 2 without measurement errors (σ=0).
Fig. 6. Comparison between the desired and estimated heating power for different measurement errors in example 2 when r=24.
Table 1. Physical and thermophysical properties (Liu and Cheng, 2008). Epidermis
Dermis
Subcutaneous
Blood
c (J/kg K)
3600
3400
3060
3770
ωb (kg/m3)
0.0
1.5
1.25
-
k (W/m K)
0.26
0.52
0.21
-
ρ (kg/m3)
1200
1200
1000
1060
Table 2. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 1. r=8
r=12
r=16
r=20
r=24
r=28
r=32
r=36
r=40
σ=0
3.172
3.110
2.714
2.491
2.425
2.476
2.618
2.841
3.125
σ=0.1
4.783
4.277
4.025
3.831
3.524
3.241
3.005
2.746
2. 941
σ=0.2
5.898
5.439
5.129
5.009
4.808
4.433
4.102
3.781
3.644
Table 3. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 2. r=8
r=12
r=16
r=20
r=24
r=28
r=32
r=36
r=40
σ=0
2.040
1.895
1.601
1.441
1.394
1.397
1.403
1.425
1.496
σ=0.1
3.469
3.171
2.744
2.351
2.059
1.993
1.815
1.802
1.805
σ=0.2
4.112
3.625
3.191
2.948
2.764
2.460
2.385
2.377
2.506
Highlights
•
The unknown time-dependent part of heat source was estimated by using the sequential method from the temperature measurements at one point.
•
One dimensional Pennes equation was used for modeling the heat transfer in a multilayer tissue, which is stratified into epidermis, dermis, and subcutaneous layers.
•
The direct problem was solved by using the finite difference method.
•
Two test cases were presented to examine the ability of the inverse analysis on the estimation.
•
The effect of measurement errors on the accuracy of the inverse analysis was investigated.
www.elsevier.com/locate/jtherbio
PII: DOI: Reference:
S0306-4565(15)30065-6 http://dx.doi.org/10.1016/j.jtherbio.2015.07.006 TB1659
To appear in: Journal of Thermal Biology Received date: 2 March 2015 Revised date: 14 July 2015 Accepted date: 14 July 2015 Cite this article as: M. Baghban and M.B. Ayani, Source term prediction in a multilayer tissue during hyperthermia, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2015.07.006 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.
Source term prediction in a multilayer tissue during hyperthermia
M. Baghban, M. B. Ayani Department of Mechanical Engineering Faculty of Engineering Ferdowsi University of Mashhad Mashhad, Iran P.O. Box 91775-1111 Phone No.: 0098 9155117671 Fax: 0098 513 8763304 Email: [email protected]
Abstract: One of the major challenges in the use of hyperthermia to treat cancer is determining the desired heating power of external source in such a way that the thermal injury is confined to the unhealthy tissue. In this study, an inverse method based on the sequential method is proposed to estimate the desired heating power as a function of time for a successful hyperthermia treatment. In order to simulate the measured temperature, the direct problem is solved for a multilayer skin tissue to obtain the temperature data at the skin surface. These data are employed in the inverse problem to estimate the heating power of external source. Two examples are considered to examine the accuracy of the inverse analysis. In addition, the effect of measurement errors is investigated. Results show that the proposed inverse algorithm is able to determine the desired heating power of external source accurately, even in the presence of measurement errors. However, for noisy data, more temperature measurements are required to achieve reliable results.
Keywords: Inverse problem, Sequential method, Hyperthermia, Multilayer skin tissue, Heating power, Heat source.
1. Introduction Hyperthermia is a method of cancer treatment in which the temperature of tumors increased in the range of 42°C–46 °C for a specified period of time to destroy them while the healthy parts are not destroyed (Field and Bleehen, 1979; Overgaard, 1989). In the clinical application, hyperthermia can be used in two different ways, local hyperthermia and regional hyperthermia. In the local hyperthermia, very high temperatures are applied to destroy the tumor while in the regional hyperthermia the temperatures increase a few degrees to help other cancer treatments such as chemotherapy and radiotherapy (Hildebrandt et al., 2002).
Nomenclature
∆x Y Z
specific heat, J/kg K objective function heat transfer coefficient, W/m2 oC thermal conductivity, W/m oC length, m number of nodes at spatial coordinate number of temporal measurements number of future time steps source term, W/m3 time-dependent part of heat source, heating power, W/m2 temperature, K, oC temporal coordinate, s increment of time domain, s spatial coordinate, m increment of spatial coordinate, m measured temperature, K, oC sensitivity matrix
Greek symbols β
scattering coefficient, 1/m
c F(q) h k L M N r S(x,t) q(t) T(x,t) t ∆t
x
ρ ζ σ ω
density random variable standard deviation of the measurement errors blood perfusion
Subscript b es d E f if s P W ∞
blood estimated desired east node final time interface s-th layer central node west node cooling fluid
Superscript n *
time level Initial guess
To have a successful thermal therapy, it is necessary to control the temperature. For this goal, many researchers attempted to solve inverse bio-heat conduction in estimating the boundary condition and heat source. Loulou and Scott (2002) used the conjugate gradient method to optimize the heat source that results in the desired thermal dose. Kuznetsov (2006) applied the minimum principal of Pontryagin and determined the optimal volumetric heating rate. Yang (2011) estimated the boundary condition in hyperbolic bio-heat conduction by using the sequential method. Dhar et al. (2012) investigated analytically optimal distribution of time dependent heating power in a multilayer tissue. Lee et al. (2013) developed an inverse hyperbolic analysis for estimating the unknown time-dependent surface heat flux of a living skin tissue. Yang (2014a) predicted the boundary condition in two-dimensional non-Fourier multilayer bioheat conduction. Lee et al (2014) used the conjugate gradient method to estimate the unknown time-dependent surface heat flux in a multilayer tissue based on the hyperbolic model of bio-heat
conduction. Jalali et al. (2014) used the conjugate gradient method to estimate simultaneously the heat source and heat transfer coefficient at the skin surface for the desired temperature distribution. Yang (2014b) determined the heat strength that destroys tumor tissue but does not damage normal tissue. This study deals with determining the unknown time-dependent heating power of external source using the sequential algorithm of the inverse bio-heat transfer. In order to simulate the measured temperature, the direct problem is solved for a one-dimensional multilayer tissue to obtain the temperature data at the skin surface. The accuracy of the inverse method is checked by utilizing two examples.
2. Direct problem The physical model applied in this work is sketched in Fig. 1. Assume that the thermal properties are constant and the metabolic heat generation is negligible compared to the heat source due to external heating; the thermal model of bio-heat transfer in the s-th layer of the tissue is given as (Pennes, 1948): ρscs
∂T ∂t
= ks
∂2T ∂x 2
+ ωbscbs (Tb −T ) + S (x ,t )
(1)
where ρ, c and k are the density, specific heat and thermal conductivity of the tissues, respectively. ωb, cb and Tb denote the blood perfusion rate, specific heat and the arterial blood temperature, respectively. S(x,t) is the heat source term where heat can be generated from external sources such as laser or ultrasound. The initial condition and boundary conditions are expressed as:
T (x ,t ) = 37 °C for t = 0
(2)
∂T (x , t ) = h (T ∞ −T (x , t ) ) ∂x
−k
(3)
at x = 0
T (x , t ) = 37 °C at x = L
(4)
where T∞ and h represent the temperature and heat-transfer coefficient of the cooling fluid on the skin surface, respectively. At the interface of two layers, the temperature and heat flux are continuous. The boundary conditions at the interface of i-th and j-th layers (x=xif) are as follows:
T i (x , t ) = T j (x , t ) −k i
(5)
at x = x if
∂T j ( x , t ) ∂T i (x , t ) = −k j ∂x ∂x
(6)
at x = x if
The fully implicit finite difference method (Ferzige and Peric, 2002) is used to solve the direct problem, subject to the initial condition and the boundary conditions. Therefore, the governing equation (Eq. (1)) can be expressed as the following: k s − (∆x )2
ρ c 2k s T E n +1 + s s + + ωbs cbs ∆t (∆x )2
ρs c s n T P ∆t
k s T P n +1 − (∆x )2
TW n +1 =
(7)
+ (ωbs cbs )Tb + S (x ,t )
where TPn+1 is the temperature of node P at the time n+1 (see Fig. 2). ∆x and ∆t are the increments of spatial coordinate and time domain, respectively. TE and TW are the temperature of the east and west nodes, respectively. In the direct problem, the source term is known, while in the inverse problem, it is considered as an unknown value and needs to be predicted by measuring the temperatures at the skin surface.
3. Inverse problem
In this study, the sequential method combined with the concept of the future time (Beck et al., 1985) is used to estimate the heat source term. According to the well-known Beer's Law, the heat source term is given by (Deng and Liu, 2002),
S (x ,t ) = β q (t )e −βx
(8)
where β is the scattering coefficient. q(t) is the time-dependent part of source term, which is the heating power of the external source. The objective of the inverse analysis is to estimate q(t) from the knowledge of temperature measurement taken at the skin surface (x=0). To estimate the transient heating power, the heating power is broken up into ‘‘n’’ elements. Therefore, the governing equation at t=tm is rewritten by the following equations: ρscs
∂T (x ,t m ) ∂t
= ks
∂2T (x ,t m ) ∂ 2x
+ ωbscbs (Tb −T (x ,t m )) + qm +1β e −βx
(10)
T (x , t ) = T m −1 for t = t m −1
∂T (x , t m ) = h (T ∞ −T (x , t m ) ) ∂x
−k
T (x , t m ) = 37 °C
(9)
at x = 0
(11)
(12)
at x = L
At the interface of i-th and j-th layers (x=xif), the boundary conditions are as follows:
T i (x , t m ) = T j (x , t m ) −k i
(13)
at x = x if
∂T j (x , t m ) ∂T i ( x , t m ) = −k j ∂x ∂x
at x = x if
(14)
The temperature distribution at t=tm is determined by solving Eqs. (9)-(14). In order to estimate qm+1, it is supposed that the values q1, q2, ..., qm have been previously calculated. In addition, the heating power components at the future time are assumed to be constant over "r" future time steps, i. e,
(15)
qm+2 = ... = qm+r −1 = qm +r = qm +1
The objective function for the identification of qm+1 is expressed by the sum of square residuals between the measured and estimated temperatures as follows: r
F (qm+1) = ∑[Y m+i −Tm+i (qm+1)]
2
(16)
i =1
Equation (16) is minimized with respect to the qm+1 to obtain, r dF = −2∑ [Y m+i −Tm+i (qm+1)][ Z i ] dqm+1 i =1
(17)
where the step sensitivity coefficients, Z, at the sensor location defined as:
∂T ∂Tm+2 ∂Tm+r Z = m+1 … ∂qm+1 ∂qm+1 ∂qm+1
(18)
Employing a Taylor series expansion of Tm+1(qm+1) about an initial guess, q*, qm+1 is calculated as: (19)
r
∑(Y qm +1 = q* + i =1
m +i
−Tm +i )Z i
r
∑Z
2 i
i =1
In linear problems such as the present study, Eq. (19) is an accurate equation for estimating the unknown heating power.
3.1 Sensitivity Problem To calculate the sensitivity coefficients (Eq. (18)), the direct problem given by Eqs. (9)-(14) is differentiated with respect to qm+1 as: ρscs
∂Z (x ,t m ) ∂t
Z (x , t m ) = 0
= ks
∂2Z (x ,t m )
for
∂2x
−ωbscbs Z (x ,t m ) + β e −βx
t m = t m −1
(20)
(21)
k
∂Z (x , t m ) = h Z (x , t m ) at x = 0 ∂x
(22)
Z (x , t m ) = 0 at x = L
(23)
Z i (x , t m ) = Z j (x , t m ) at x = x if
(24)
−k i
∂Z j ( x , t m ) ∂Z i (x , t m ) = −k j ∂x ∂x
at x = x if
(25)
The finite difference method is used to solve the sensitivity problem.
4. Results and discussion The sequential algorithm developed in the present study is tested by using it for the problem of estimating the time-dependent heating power, q(t), in a one-dimensional multilayer tissue. Since experimental data are not available, a real heat source is assumed to calculate the temperatures at the measurement location (x=0). The results are considered as the computed temperature, Y. In order to simulate the measured temperature containing measurement errors, random errors of standard deviation, σ, are added to the temperature computed from the solution of the direct problem. It can be shown in the following equation,
Y =Y ex +σζ
(26)
where ζ is a random variable within -2.576 to 2.576 for a 99% confidence bounds. The physical and thermophysical properties of all layers are listed in Table 1 (Liu and Cheng, 2008). The skin surface is a convection boundary with the convection coefficient h=200 W/m2 oC and the cooling temperature T∞=5 oC. The scattering coefficient for all layers is taken as β=100 m-1 (Dhar et. al, 2012). The final time is chosen as tf =1000s. The direct solution for a one-layer tissue is validated by Deng's analytical solution (Deng and Liu, 2002) by assuming the same conditions. To demonstrate the accuracy and efficiency of the inverse solution in predicting the
time dependent heating power, two different functions involving abrupt changes are considered in this study. In all examples, the estimated time-dependent heating power, qes(t), obtained with an initial guess q*(t)=0.
Example 1. The desiered transient heating power, qd (t), is assumed in the following form,
500 qd (t ) = 1000 500
tf 3 tf 2t ≤t ≤ f 3 3 2t f < t ≤ tf 3 0 ≤t <
(27)
Fig. 3 displays the effect of future time steps number on the accuracy of estimation without measurement errors. It can be seen clearly that increasing the number of future time steps, improves the estimating stability but decreases the accuracy. To determine the accuracy of the inverse solution, the relative percentage error is defined by: N
qes −qd ×100, qd m=1
ERR = ∑
(28)
where N is the number of temporal measurements. The relative percentage errors for different number of future time steps and measurement errors are listed in Table 2. It is obvious from the results that r=24 could present a stable and accurate estimation when σ=0. In addition, it can be concluded that the heat source term can be estimated accurately, even with noisy data. However, for noisy data, more number of future time steps is required to get an accurate solution. The
influence of the measurement errors on the solution for example 1, when the future time steps number is 24, is shown in Fig. 4. For r=24, the relative errors for q(t) are computed as 2.425%, 3.524% and 4.808%, for σ=0, σ=0.1 and σ=0.2,
respectively. The results display that the
proposed solution is an accurate and stable method to determine the time-dependent heating power.
Example 2. The desired time-dependent heating power, qd (t), is taken in the following form,
t 0 ≤t ≤ f 500 4 2000(t − t f )/ t + 500 t f ≤ t ≤ t f f 4 4 2 qd (t ) = −2000(t − 3t f )/ t + 500 t f ≤ t ≤ 3tf f 4 2 4 3t f 500 ≤t ≤ tf 4
(29)
The effects of future time steps and measurement errors on the accuracy and stability of estimation is demonstrated in Figs. 5 and 6, respectively. In addition, the relative percentage errors for different number of future time steps and measurement errors are listed in Table 3. It can be seen, the relative percentage error is decreases to the minimum value and then increases when the future time steps number increases. However, increasing r results in a more stable solution. Also, it is clear that the relative error is increased with the increase of the measurement errors. For r=24, the relative error decreases 49% as the measured error reduces from σ=0.2 to
σ=0. The results show that the estimated values of the time-dependent heating power agree well with the desired values even for the case with the measurement errors.
5. Conclusion
An inverse algorithm based on the sequential method was used to estimate the time-dependent part of source term in a one-dimensional multilayer skin tissue from the temperature measurement at the skin surface. The Pennes bio-heat equation was applied to simulate thermal behavior in the tissue. Two test cases were presented to check the ability and accuracy of the proposed inverse solution. The effect measurement noise on the accuracy of the estimation was studied. The results showed that the agreement between the estimated and desired values of the heat source was excellent, even with noisy data.
References Beck, J. V., Blackwell, B., Clair Jr, C. R. S., 1985. Inverse heat conduction: Ill-posed problems. Wiley, New York. Deng, Z. S., Liu, J., 2002. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. Journal of biomechanical engineering, 124 (6), 638-649. Dhar, P., Dhar, R., Dhar, R., 2012. Analytical study on optimal distribution of heating power in hyperthermia. International Communications in Heat and Mass Transfer, 39 (3), 419-423. Ferziger, J. H., Peric¸ M., 2002. Computational Methods for Fluid Dynamics. Springer, New York. Field, S. B., Bleehen, N. M., 1979. Hyperthermia in the treatment of cancer. Cancer treatment reviews, 6 (2), 63-94. Wust, P., Hildebrandt, B., Sreenivasa, G., Rau, B., Gellermann, J., Riess, H., Felix, R., Schlag, P. M., 2002. Hyperthermia in combined treatment of cancer. The lancet oncology, 3 (8), 487-497.
Jalali, A., Ayani, M. B., Baghban, M., 2014. Simultaneous estimation of controllable parameters in a living tissue during thermal therapy. Journal of thermal biology, 45, 37-42. Kuznetsov, A. V., 2006. Optimization problems for bioheat equation. International Communications in Heat and Mass Transfer, 33 (5), 537-543. 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), 2630-2643.
Lee, H. L., Chen, W. L., Chang, W. J., Yang, Y. C., 2015. Estimation of surface heat flux and temperature distributions in a multilayer tissue based on the hyperbolic model of heat conduction. Computer methods in biomechanics and biomedical engineering, 18 (14), 15251534. Liu, K. C., Cheng, P. J., 2008. Finite propagation of heat transfer in a multilayer tissue. Journal of Thermophysics and Heat Transfer, 22 (4), 775-782. Loulou, T., Scott, E. P., 2002. Thermal dose optimization in hyperthermia treatments by using the conjugate gradient method. Numerical Heat Transfer: Part A: Applications, 42 (7), 661-683. Overgaard, J., 1989. The current and potential role of hyperthermia in radiotherapy. nternational Journal of Radiation Oncology • Biology • Physics, 16 (3), 535-549. Pennes, H. H., 1948. Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of applied physiology, 1 (2), 93-122. Yang, C. Y., 2011. Boundary estimation of hyperbolic bio-heat conduction. International Journal of Heat and Mass Transfer, 54 (11), 2506-2513.
Yang, C. Y., 2014a. Boundary prediction of bio-heat conduction in a two-dimensional multilayer tissue. International Journal of Heat and Mass Transfer, 78, 232-239. Yang,
C.
Y.,
2014b.
Determining
the
heat
strength
required
in
hyperthermia
treatments. International Communications in Heat and Mass Transfer, 57, 282-285.
Figure Captions: Fig. 1. Schematic of a one-dimensional multilayer tissue. Fig. 2. Schematic of a node. Fig. 3. Comparison between the desired and estimated heating power for different future time steps in example 1 without measurement errors (σ=0). Fig. 4. Comparison between the desired and estimated heating power for different measurement errors in example 1 when r=24. Fig. 5. Comparison between the desired and estimated heating power for different future time steps in example 2 without measurement errors (σ=0). Fig. 6. Comparison between the desired and estimated heating power for different measurement errors in example 2 when r=24.
Table Captions: Table 1. Physical and thermophysical properties (Liu and Cheng, 2008). Table 2. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 1. Table 3. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 2.
Fig. 1. Schematic of a one-dimensional multilayer tissue.
Fig. 2. Schematic of a node.
Fig. 3. Comparison between the desired and estimated heating power for different future time steps in example 1 without measurement errors (σ=0).
Fig. 4. Comparison between the desired and estimated heating power for different measurement errors in example 1 when r=24.
Fig. 5. Comparison between the desired and estimated heating power for different future time steps in example 2 without measurement errors (σ=0).
Fig. 6. Comparison between the desired and estimated heating power for different measurement errors in example 2 when r=24.
Table 1. Physical and thermophysical properties (Liu and Cheng, 2008). Epidermis
Dermis
Subcutaneous
Blood
c (J/kg K)
3600
3400
3060
3770
ωb (kg/m3)
0.0
1.5
1.25
-
k (W/m K)
0.26
0.52
0.21
-
ρ (kg/m3)
1200
1200
1000
1060
Table 2. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 1. r=8
r=12
r=16
r=20
r=24
r=28
r=32
r=36
r=40
σ=0
3.172
3.110
2.714
2.491
2.425
2.476
2.618
2.841
3.125
σ=0.1
4.783
4.277
4.025
3.831
3.524
3.241
3.005
2.746
2. 941
σ=0.2
5.898
5.439
5.129
5.009
4.808
4.433
4.102
3.781
3.644
Table 3. The relative percentage errors (ERR) with the variation of the measurement errors and the future time steps number in example 2. r=8
r=12
r=16
r=20
r=24
r=28
r=32
r=36
r=40
σ=0
2.040
1.895
1.601
1.441
1.394
1.397
1.403
1.425
1.496
σ=0.1
3.469
3.171
2.744
2.351
2.059
1.993
1.815
1.802
1.805
σ=0.2
4.112
3.625
3.191
2.948
2.764
2.460
2.385
2.377
2.506
Highlights
•
The unknown time-dependent part of heat source was estimated by using the sequential method from the temperature measurements at one point.
•
One dimensional Pennes equation was used for modeling the heat transfer in a multilayer tissue, which is stratified into epidermis, dermis, and subcutaneous layers.
•
The direct problem was solved by using the finite difference method.
•
Two test cases were presented to examine the ability of the inverse analysis on the estimation.
•
The effect of measurement errors on the accuracy of the inverse analysis was investigated.