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

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.


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][16][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.

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, 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  cl are the thickness and thermal conductivity of the casing layer, 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 −4 ∼10 −5 ).The computational algorithm for the foregoing procedure is given by the following.
(1) Give an initial guess for .

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).

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 −2 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 −4 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.

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 (0.007,  ) = 43 ∘ C and two others in the normal tissue with (0.002,  ) and (0.021,  ) = 39 ∘ 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.

Results and Discussions
In the present work, FDTD method is used for numerical solution of the direct, sensitivity, and adjoint problems.A Moreover, the computational precision () is set at 10 −5 .In the following simulation, the desired temperature distribution is determined by precisely prescribing a known .
Figure 2 indicates the sensitivity of the temperaturetime 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/m 2∘ 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, | retrieved −  prescribed |/ prescribed .
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  = 200 W/m 2∘ C as well as the desired one, which is obtained from the application of a known  at  = 1000 s.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.
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.
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.

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.

Figure 1 :
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.

( 6 )( 7 )( 8 )
Compute the conjugate coefficient   and the direction of descent   .Solve the sensitivity problem.Compute the search step size   .

Figure 2 :
Figure 2: Temperature sensitivity resulting from variations in .

Figure 3 :
Figure 3: Convergence curve for different initial guesses of .

Table 1 :
Biophysical characteristics of model.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.