Journal of Applied Mathematics

Volume 2013 (2013), Article ID 734020, 9 pages

http://dx.doi.org/10.1155/2013/734020

Research Article

An Inverse Problem of Temperature Optimization in Hyperthermia by Controlling the Overall Heat Transfer Coefficient

Seyed Ali Aghayan,¹ Dariush Sardari,¹ Seyed Rabii Mahdi Mahdavi,² and Mohammad Hasan Zahmatkesh³

¹Faculty of Engineering, Science and Research Branch, Islamic Azad University, P.O. Box 14515-775, Tehran, Iran

²Medical Physics, Tehran University of Medical Science, Tehran, Iran

³Medical Physics, Shahid Beheshti University, Tehran, Iran

Received 14 January 2013; Revised 1 July 2013; Accepted 8 July 2013

Academic Editor: Tin-Tai Chow

Copyright © 2013 Seyed Ali Aghayan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A novel scheme to obtain the optimum tissue heating condition during hyperthermia treatment is proposed. To do this, the effect of the controllable overall heat transfer coefficient of the cooling system is investigated. An inverse problem by a conjugated gradient with adjoint equation is used in our model. We apply the finite difference time domain method to numerically solve the tissue temperature distribution using Pennes bioheat transfer equation. In order to provide a quantitative measurement of errors, convergence history of the method and root mean square of errors are also calculated. The effects of heat convection coefficient of water and thermal conductivity of casing layer on the control parameter are also discussed separately.

1. Introduction

Hyperthermia treatment is recognized as an adjunct cancer therapy technique following surgery, chemotherapy, and radiation techniques. In hyperthermia, the tumor cells will be overheated, typically 40–45°C, to damage or kill the cancer cells [1, 2]. The determination of temperature distribution throughout the biological tissue by solving the simple Pennes bioheat Equation (Eq) and increment of tumor temperature to therapeutic value avoiding the overheating of the healthy tissue are important issues of investigation in hyperthermia [3]. Because of obvious reasons, which are not going to be discussed here, the water (or any other dielectric liquid) bolus is known as a necessary element of all contact applicators used in clinics in order to provide local heating of malignant tumors by means of electromagnetic (EM) energy [4]. There are numerous papers discussing the effect of the water coupling bolus on the temperature of surface and tissue specific absorption rate (SAR) during hyperthermia treatment [5, 6]. Sherar et al. [7] presented a new way to design the water bolus for external microwave applicators which increased the available treatment field size by removing the central hotspot caused by the increased power from the applicator in this region. They also used beam shaping bolus with an array of saline-filled patches inside the coupling water bolus for superficial microwave hyperthermia [8]. Ebrahimi-Ganjeh and Attari [9] numerically computed and compared both the SAR distributions and effective field size in the presence and absence of water bolus. Neuman et al. [10] examined the effect of the various thickness of water bolus coupling layers on the SAR patterns from dual concentric conductor based on the conformal microwave array superficial hyperthermia applicators.

One aspect of our investigation is optimal control theory, which is defined as a mathematical approach of manipulating the parameters affecting the behaviour of a system to produce the desired or optimal outcome. A recent development in this field is the optimal control of systems with distributed parameters which especially include the biological tissue in course of the hyperthermia treatment planning [11]. Butkovasky [11] treated the fundamentals in the theory of optimal control of systems with distributed parameters. Wagter [12] studied an optimization procedure to calculate the transient temperature profiles in a plane tissue by multiple EM applicators. Kuznetsov [13] used the minimum principle of Pontryagin to solve an optimum control problem for heating a layer of tissue in order to destroy a cancerous tumor. Although the Levenberg-Marquardt and simplex methods are common techniques for solving optimization problems, we tend to develop a methodology to numerically optimize hyperthermia treatments based on an inverse heat conduction problem by using the conjugate gradient method with an adjoint equation (CGMAE). One of the advantages of the CGMAE is its insensitivity to initial variables; however, it requires a numerical method for each iteration to determine the temperature fields, which involves computing the gradients of objective function. In the last two decades, CGMAE has been widely used in the solution of the general inverse heat transfer problems (IHTPs) regarding the optimization procedure [14] especially in therapeutic applications [15–17]. The most common approach for determining the optimal estimation of the model parameters, functions, or boundary conditions is to minimize an objective function.

In this study, the objective function to be minimized is based on the criterion of minimal square root of the difference between the desired temperature and the computed one. We demonstrate that optimum heating conditions inside the biological tissue, which is caused by an applied radiofrequency (RF) heat source, can be numerically attained within the total time of the process by means of controlling the overall heat transfer coefficient associated with the cooling system. To do this and to reach the appropriate configuration of cooling system, thermal properties of water bolus and casing layer are considered as variable parameters. Recently, few studies have been conducted with the aim of beam-forming in hyperthermia which are based on utilizing the saline patches inserted in the cooling system or the solid absorbing blocks between the applicator and water bolus [7, 8, 18]. However, best to our knowledge, none of them has used an inverse model to control the thermal or electrical properties of these additional patches or blocks. In the perspective of the model, we believe that the desired temperature profile inside the tumor according to its shape and size can be achieved by inversely optimizing the shape, thickness, electrical, and thermal properties of solid or liquid patches which are inserted in the cooling system. Therefore, the proposed algorithm developed in this study could be effectively used by physicians to design an optimal treatment plan in a large variety of therapeutic treatments.

2. Optimization Scheme

2.1. Governing Equations and Model Geometry

Schematic diagram of the studied control zone is shown in Figure 1. In order to simplify demonstrating the methodology, heat transfer model is assumed to be one-dimensional. It should be noted that although we could develop a model with an irregular shape of a tumor embedded in the healthy tissue, a homogeneous slab of tissue, including the tumor region, is considered. Temperature distribution inside the tissue could be obtained by solving the Pennes bioheat equation [19] as follows:

where is the overall heat transfer coefficient, which defines the effectiveness of the cooling system. Because the temperature distribution is very sensitive to the heat transfer between the cooling liquid and the so-called tissue, its numerical value must be properly estimated. It is shown that the expression of is [20] as follows:

Here, is the heat convection coefficient of water and and are the thickness and thermal conductivity of the casing layer, respectively.

Figure 1: Schematic diagram of control zone. Cooling system consists of water bolus embedded in a casing layer which is placed between the applicator and biomedia. Regions I and II represent the healthy tissue and tumor region, respectively.

For the inverse problem proposed here, the overall heat transfer coefficient of cooling system is considered as an unknown function, while the other quantities within the formulation of the direct problem presented previously are assumed to be known precisely. Moreover, the temperature distribution inside the tissue region is considered available.

The optimization procedure includes maximizing the temperature inside the tumor while minimizing it in the normal surrounding tissue via a proposed objective function. In other words, we tend to achieve the desired temperature,, at points whether in the tumor region or in the healthy one, by optimal control of within the total time of the treatment process. Therefore, the following function must be minimized:

where is the computed temperature. In order to minimize under the constraints mentioned in (1)–(4), we introduce a new function called the “Lagrangian function” which associates (6) and (1) with the following equation:

where is the Lagrange multiplier obtained by solving the adjoint problem which is defined by the following set of equations:

where is the Dirac function. Adjoint problem can also be solved by the same method as the direct problem of (1), that is, by using the finite difference time domain method (FDTD). It can be shown that the gradient of residual functional, , could be expressed by the following analytical expression:

The iterative procedure of CGMAE for the optimization of is given by

Here, is the step size and is the descent direction of the th iteration which is defined as follows:

The coefficient of conjugate gradient, , is given by

Then, the step size could be derived from the following equation:

Here, is the answer to the sensitivity problem. As mentioned in the appendix, we present the methodology used to analytically derive the sensitivity, the adjoint problems, and the gradient equation. The iteration procedure is repeated until it satisfies the following truncation criteria and then leads to optimal control of the overall heat transfer coefficient as follows:

where is a small number (10⁻⁴~10⁻⁵). The computational algorithm for the foregoing procedure is given by the following.

(1)Give an initial guess for .(2)Solve the direct problem.(3)Examine (14), and if convergence is satisfied, stop; otherwise, continue.(4)Solve the adjoint problem.(5)Compute the gradient equation.(6)Compute the conjugate coefficient and the direction of descent .(7)Solve the sensitivity problem.(8)Compute the search step size .(9)Increment the controlling parameter using (10) and return to Step 1.

2.2. Model Parameters

2.2.1. Electromagnetic Model

The spatial and temporal behaviors of the EM field inside the tissues and are governed by Maxwell’s equations. In a homogeneous tissue, these equations are expressed by [21] the following:

where and are, respectively, complex electric and magnetic fields, is the angular frequency, and , and are the electrical permittivity, the electrical conductivity, and the magnetic permeability of tissue, respectively. is the complex permittivity which is defined by

The volumetric EM power deposition, which averaged over time, is calculated by

where is the root mean square magnitude [V/m] of the incident electric field in a point which could be calculated through (15).

2.2.2. Cooling System and Tissue Properties

Thermophysical properties of tissue and blood, pertaining to the model for numerical simulation, are provided in Table 1 [3, 22, 23]. The thickness of water is set at 10⁻² m and the electrical permittivity and electrical conductivity of water, related to an RF applicator operating at 432 MHz, are 78 and 0.0300 S/m, respectively [23]. Moreover, thicknesses of 0.03 m and 1.5 10⁻⁴ m are used, respectively, for tissue and casing layer which are typically made of thin flexible Silicon layers. Temperature of the cooling system, , is set at 15°C.

Table 1: Biophysical characteristics of model.

Component
Parameter
Value
Units

Blood

4200
j/kg °C

1000
kg/m³

0.0005
mL/s/mL

Tissue

3770
j/kg °C

0.5
W/m °C

33800
W/m³

37
°C

37
°C

34
°C

55
—

998
kg/m³

1.1
S/m

3. Assumptions

In order to simplify demonstrating the methodology, heat transfer is assumed to be one-dimensional with a constant rate of tissue metabolism and blood perfusion. The target temperatures, after the elapse of allowed heating period, are set between 42 and 45°C for tumoral region and 42°C or lower for normal tissues. In addition, we consider three different nodes: one in the tumor region with °C and two others in the normal tissue with and °C as the desired temperatures. As the thickness of the casing layer is much smaller than its wavelength, it could be neglected in the EM analysis. Moreover, the thermal contact between the casing layer and tissue is assumed to be perfect so that the thermal contact resistance between these layers is negligible.

4. Results and Discussions

In the present work, FDTD method is used for numerical solution of the direct, sensitivity, and adjoint problems. A uniform grid of 31 × 31 mesh in space and time is utilized for all the results presented hereafter. The optimization algorithm presented previously is programmed in Fortran 90, compiled through the software Compaq Visual Fortran Professional Edition on a computer with an Intel Core i7 3.5 GHz and 16 Gb of RAM. The final time,, is chosen equal to 1800 s. Moreover, the computational precision () is set at 10⁻⁵. In the following simulation, the desired temperature distribution is determined by precisely prescribing a known .

Figure 2 indicates the sensitivity of the temperature-time curve for variations in for three different values of . Because of the large values of the sensitivity with respect to , obviously the proposed optimization model is “well-conditioned.” To illustrate the accuracy of the present inverse analysis, the initial approximations of were either increased or decreased from its known value. Figure 3 shows the number of convergent iterations with different initial guesses for . Changing the initial value, for example, by 15% and 30% of its known value does not affect the solution and the optimization object is still satisfied. Moreover, it is clear that the number of convergent iterations is reduced when the initial guess value varies within a 15% deviation from the known one. In addition, this figure shows that the minimization of the objective function is valid when varies between 125 and 200 W/m2°C. In order to provide a quantitative measurement of errors, the convergence history of the presented model is also calculated. Figure 4 compares the convergence history is of three different initial guesses. Convergence history is measured by the relative error at each iteration, that is, .

Figure 2: Temperature sensitivity resulting from variations in .

Figure 3: Convergence curve for different initial guesses of .

Figure 4: Comparison of convergence histories of CGM for different initial guesses— : 125, □: 150, and ○: 200 W/m2°C for .

The computed temperature distribution is obtained when the minimization of the objective function is satisfied, that is, when the optimal value of is achieved. Figure 5 shows the computed temperature distribution with initial guess of  W/m2°C as well as the desired one, which is obtained from the application of a known at . It is observed that the recovered temperature distribution is in good agreement with the desired one. In addition, the computed RMS error is equal to 0.24. Figure 6 presents the desired temperature distributions result from the solution of the nonstationary thermal problem during the total heating process. It is found to be sufficient to elevate the temperature inside the tumor to 42°C within approximately 20 minutes, while the temperature on the surface remains approximately constant in this period of time.

Figure 5: A comparison between recovered and desired temperature distribution at  s.

Figure 6: Transient spatial temperature distribution inside the tissue.

According to the foregoing discussion, as (5), and in order to enhance the performance of the cooling system, the thermal conductivity of casing layer and the heat convection coefficient of water have been considered as variable parameters. Although the casing layer is very thin, as shown in Figure 7, the influence of thermal conductivity of casing layer is remarkable. It could be observed that a slight increase in thermal conductivity reduces significantly the temperature at the skin surface. In addition, the position of the maximum temperature is shifted slightly toward the skin. The second possibility to control the effectiveness of the cooling system is through the variation of the heat convection coefficient of water. Figure 8 indicates the influence of two different values of on the temperature of tissue. It could be observed that a 75% increase in results in a decrease of about 2.5°C in the skin and only 2°C in the maximum. The same results could be obtained by only 20% of the variation in the thermal conductivity of casing layer. However, depending on the technical facilities, one can make an optimal decision in this case.

Figure 7: Temperature distribution at  s for two different thermal conductivities of casing layer. Solid line:  W/m°C,  W/m2°C; dashed line:  W/m°C,  W/m2°C.

Figure 8: Temperature distribution at  s for two different heat convection coefficients of water. Solid line:  W/m2°C, W/m2°C; dashed line:  W/m2°C,  W/m2°C.

Although the thickness and the electrical conductivity of water can also be optimized in the treatment process, we did not consider them in this study and we mainly focused on the influence of the overall heat transfer coefficient on the temperature field inside the tissue. In the perspective of the present model, we believe that the desired temperature profile inside the tumor, according to its shape and size, can be achieved by inversely optimizing the shape, thickness, electrical, and thermal properties of solid or liquid patches which are inserted in the cooling system. Such an optimization problem should also incorporate the best configuration of the additional patches which reasonably must resemble the tumor shape, for example, L-form patches for L-form tumors. However, the clinician will first scan the geometry of the tumor and measure the tumor size and depth using an MRI or CT or other noninvasive imaging techniques. Although the present model is based on the 1D heat transfer, but with some modifications in the governing formulations, such as , it has the potential to provide a more comprehensive understanding of the cooling system ability for the purposes of beam-forming in two or three dimensions. This modified optimization algorithm is then implemented to the tumor and its surrounding tissue to determine the optimal properties of cooling system.

5. Conclusions

In this paper, we demonstrated the validity of the CGMAE in the optimization of the temperature distribution inside the tissue by optimal control of the overall heat transfer coefficient associated with the cooling system. The temperature distribution inside the tissue, based on the bioheat transfer equation, has also been numerically computed. In order to evaluate the effectiveness of the cooling system, the thermal conductivity of casing layer and the heat convection coefficient of water have been considered separately. According to results, the presented model allows obtaining the best properties of water and casing layer in the light of optimization and therefore significantly improves the treatment outcome. However, this method shows great promise for further extensions in both model and clinical performances for any shaped targeted tumor in hyperthermia treatments.

Appendix

In order to find the minimum of the residual function by gradienttype method, it is required to derive its gradient by the following processes first, computation of the variation problem when the unknown has a small variation, and second, inspection of the conditions under which we obtain a stationary point for the Lagrange functional. To compute the variation problem, the unknown function is perturbed by ; that is,

Then, the temperature undergoes the variation as follows:

where refers to the perturbed variables. By inserting the a aforementioned equations in the direct problem and subtracting the original (1)–(4) from the resulting expressions and ignoring the second order terms, the following direct problem in variation is obtained:

The second step is to calculate the variation of the objective function, , which is given by

We obtain the variation of the augmented functional, , by considering the variation problem given by and the variation of the objective function, , as follows:

Using integration by parts in ,

and by using the boundary conditions of the direct problem in variations (Equations – and the rearrangement of the mathematical terms) the final form of can be written as follows:

To satisfy the optimality condition, that is, , it is required that all coefficients of are to be vanished. Finally, the following term is left:

Since we know that and are, respectively, the direct and the adjoint problem solutions, it follows that

and then

Based on the definition, the variation of the functional can be written as follows:

where the associated inner scalar product of two functions and in the Hilbert space of is defined by

From the definition of the scalar product, the variation of the functional can be written as follows:

A comparison between , , and reveals that the gradient of the residual functional is given by (9). More details about the derivation of the equations are mentioned in [24–26].

Nomenclature

| Heat capacity
: | Vector of descent direction
: | Computational precision
: | Complex electric field
: | Residual functional
: | Residual functional variation
: | Gradient of residual functional
: | Heat convection coefficient
: | Complex magnetic field
: | Thermal conductivity
: | Thickness of casing layer
: | Lagrangian functional
: | Number of nodes
: | Tissue metabolism
: | Volumetric heat source
: | Temperature
: | Time
: | Overall heat transfer coefficient
: | Space
: | Desired temperature.

Greek Symbols

: | Descent parameter
: | Descent direction parameter
: | Dirac function
: | Electrical permeability
: | Density
: | Magnetic permeability
: | Electrical conductivity
: | Lagrange multiplier
: | Angular frequency
: | Blood perfusion.

Subscripts

: | Arterial
: | Blood
: | Core
: | Casing layer
: | Final
: | Initial time
: | Tissue
: | Ambient.

Superscripts

:
Iteration number.

References

S. R. H. Davidson and D. F. James, “Drilling in bone: modeling heat generation and temperature distribution,” Journal of Biomechanical Engineering, vol. 125, no. 3, pp. 305–314, 2003.

W. G. Sawyer, M. A. Hamilton, B. J. Fregly, and S. A. Banks, “Temperature modeling in a total knee joint replacement using patient-specific kinematics,” Tribology Letters, vol. 15, no. 4, pp. 343–351, 2003.

R. Dhar, P. Dhar, and R. Dhar, “Problem on optimal distribution of induced microwave by heating probe at the tumor site in hyperthermia,” Advanced Modeling and Optimization, vol. 13, no. 1, pp. 39–48, 2011.

E. A. Gelvich and V. N. Mazokhin, “Resonance effects in applicator water boluses and their influence on SAR distribution patterns,” International Journal of Hyperthermia, vol. 16, no. 2, pp. 113–128, 2000.

M. de Bruijne, T. Samaras, J. F. Bakker, and G. C. van Rhoon, “Effects of waterbolus size, shape and configuration on the SAR distribution pattern of the Lucite cone applicator,” International Journal of Hyperthermia, vol. 22, no. 1, pp. 15–28, 2006.

E. R. Lee, D. S. Kapp, A. W. Lohrbach, and J. L. Sokol, “Influence of water bolus temperature on measured skin surface and intradermal temperatures,” International Journal of Hyperthermia, vol. 10, no. 1, pp. 59–72, 1994.

M. D. Sherar, F. F. Liu, D. J. Newcombe et al., “Beam shaping for microwave waveguide hyperthermia applicators,” International Journal of Radiation Oncology Biology Physics, vol. 25, no. 5, pp. 849–857, 1993.

J. C. Kumaradas and M. D. Sherar, “Optimization of a beam shaping bolus for superficial microwave hyperthermia waveguide applicators using a finite element method,” Physics in Medicine and Biology, vol. 48, no. 1, pp. 1–18, 2003.

M. A. Ebrahimi-Ganjeh and A. R. Attari, “Study of water bolus effect on SAR penetration depth and effective field size for local hyperthermia,” Progress in Electromagnetics Research B, vol. 4, pp. 273–283, 2008.

D. G. Neuman, P. R. Stauffer, S. Jacobsen, and F. Rossetto, “SAR pattern perturbations from resonance effects in water bolus layers used with superficial microwave hyperthermia applicators,” International Journal of Hyperthermia, vol. 18, no. 3, pp. 180–193, 2002.

A. G. Butkovasky, Distributed Control System, American Elsevier, New York, NY, USA, 1969.

C. De Wagter, “Optimization of simulated two-dimensional temperature distributions induced by multiple electromagnetic applicators,” IEEE Transactions on Microwave Theory and Techniques, vol. 34, no. 5, pp. 589–596, 1986.

A. V. Kuznetsov, “Optimization problems for bioheat equation,” International Communications in Heat and Mass Transfer, vol. 33, no. 5, pp. 537–543, 2006.

C. Huang and S. Wang, “A three-dimensional inverse heat conduction problem in estimating surface heat flux by conjugate gradient method,” International Journal of Heat and Mass Transfer, vol. 42, no. 18, pp. 3387–3403, 1999.

P. K. Dhar and D. K. Sinha, “Temperature control of tissue by transient-induced-microwaves,” International Journal of Systems Science, vol. 19, no. 10, pp. 2051–2055, 1988.

T. Loulou and E. P. Scott, “Thermal dose optimization in hyperthermia treatments by using the conjugate gradient method,” Numerical Heat Transfer A, vol. 42, no. 7, pp. 661–683, 2002.

T. Loulou and E. P. Scott, “An inverse heat conduction problem with heat flux measurements,” International Journal for Numerical Methods in Engineering, vol. 67, no. 11, pp. 1587–1616, 2006.

H. Kroeze, M. Van Vulpen, A. A. C. De Leeuw, J. B. Van De Kamer, and J. J. W. Lagendijk, “The use of absorbing structures during regional hyperthermia treatment,” International Journal of Hyperthermia, vol. 17, no. 3, pp. 240–257, 2001.

H. H. Pennes, “Analysis of tissue and arterial blood temperatures in the resting human forearm,” Journal of Applied Physiology, vol. 1, pp. 93–122, 1948.

J. P. Holman, Heat Transfer, McGraw Hill, New York, NY, USA, 9th edition, 2002.

C. Thiebaut and D. Lemonnier, “Three-dimensional modelling and optimisation of thermal fields induced in a human body during hyperthermia,” International Journal of Thermal Sciences, vol. 41, no. 6, pp. 500–508, 2002.

Z. S. Deng and J. Liu, “Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies,” Journal of Biomechanical Engineering, vol. 124, no. 6, pp. 638–649, 2002.

N. K. Uzunoglu, E. A. Angelikas, and P. A. Cosmidis, “A 432 MHz local hyperthermia system using an indirectly cooled, water-loaded waveguide applicator,” IEEE Transactions on Microwave Theory and Techniques, vol. 35, no. 2, pp. 106–111, 1987.

O. M. Alifanov, “Solution of an inverse problem of heat conduction by iteration methods,” Journal of Engineering Physics, vol. 26, no. 4, pp. 471–476, 1974.

O. M. Alifanov, Inverse Heat Transfer Problem, Springer, New York, NY, USA, 1994.

M. N. Ozisik and H. R. B. Orlande, Inverse Heat Transfer Fundamentals and Applications, Taylor & Francis, New York, NY, USA, 2000.