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

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: 𝜌_𝑡𝑐_𝑡𝜕𝑇_𝑡(𝑥,𝑡)𝑘𝜕𝑡=∇_𝑡∇𝑇_𝑡(𝑥,𝑡)+𝜔_bl𝜌_bl𝑐_bl𝑇_𝑎−𝑇_𝑡(𝑥,𝑡)+𝑄_𝑚+𝑄_𝑟(𝑥,𝑡)0<𝑥<𝑎,𝑇_𝑡(𝑥,𝑡)=𝑇_𝑐𝑥=0,−𝑘_𝑡𝜕𝑇_𝑡(𝑥,𝑡)𝑇𝜕𝑥=𝑈_∞−𝑇_𝑡(𝑥,𝑡)𝑥=𝑎,𝑇_𝑡(𝑥,𝑡)=𝑇_𝑖0<𝑥<𝑎,𝑡=0, (1)(2)(3)(4) 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: 1𝑈=1/ℎ+l/𝑘_cl. (5) Here, ℎ is the heat convection coefficient of water and 𝑙 and 𝑘_cl 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: 𝐺(𝑈)=^𝑁_𝑗=1^𝑡𝑓₀𝑇_𝑡𝑥_𝑗,𝑡,𝑈−𝑌_𝑡𝑥_𝑗,𝑡²𝑑𝑡, (6) 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: 𝐿𝑇_𝑡,Ψ_𝑡(𝑈)+,𝑈=𝐺_{𝑡𝑓0𝑎0}Ψ_𝑡(𝑥,𝑡)𝑘_𝑡𝜕²𝑇_𝑡(𝑥,𝑡)𝜕𝑥²−𝜌_𝑡𝑐_𝑡𝜕𝑇_𝑡(𝑥,𝑡)𝜕𝑡+𝜔_bl𝜌_bl𝑐_bl𝑇_𝑎−𝑇_𝑡(𝑥,𝑡)+𝑄_𝑚+𝑄_𝑟(𝑥,𝑡)𝑑𝑥𝑑𝑡, (7) where Ψ_𝑡(𝑥,𝑡) is the Lagrange multiplier obtained by solving the adjoint problem which is defined by the following set of equations: 𝑘_𝑡𝜕²Ψ_𝑡(𝑥,𝑡)𝜕𝑥²+𝜌_𝑡𝑐_𝑡𝜕Ψ_𝑡(𝑥,𝑡)𝜕𝑡−𝜔_bl𝜌_bl𝑐_blΨ_𝑡()𝑇𝑥,𝑡+2_𝑡𝑥_𝑗,𝑡−𝑌_𝑡𝑥_𝑗,𝑡×𝛿𝑥−𝑥_𝑗Ψ=00<𝑥<𝑎,_𝑡()𝑘𝑥,𝑡=0𝑥=0,_𝑡𝜕Ψ_𝑡(𝑥,𝑡)𝜕𝑥=Ψ_𝑡()Ψ𝑥,𝑡𝑈𝑥=𝑎,_𝑡(𝑥,𝑡)=00<𝑥<𝑎,𝑡=𝑡_𝑓, (8) 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: ∇𝐺(𝑈)=^𝑡𝑓₀Ψ_𝑡(𝑎,𝑡)𝑇_𝑡(𝑎,𝑡)−𝑇_∞𝑑𝑡. (9)

The iterative procedure of CGMAE for the optimization of 𝑈 is given by 𝑈^𝑙+1=𝑈^𝑙−𝛽^𝑙𝑑^𝑙𝑙=0,1,2,…. (10) Here, 𝛽^𝑙 is the step size and 𝑑^𝑙 is the descent direction of the 𝑙th iteration which is defined as follows: 𝑑^𝑙=∇𝐺^𝑙+𝛾^𝑙𝑑^𝑙−1𝑑⁰=0. (11) The coefficient of conjugate gradient, 𝛾^𝑙, is given by 𝛾^𝑙=∇𝐺𝑈^𝑙2∇𝐺𝑈^𝑙−12𝛾⁰=0. (12) Then, the step size could be derived from the following equation: 𝛽^𝑙=∑^𝑁_𝑗=1∫^𝑡𝑓₀𝑇_𝑡𝑥_𝑗,𝑡−𝑌_𝑡𝑥_𝑗,𝑡Δ𝑇_𝑡𝑥_𝑗,𝑡𝑑𝑡∑^𝑁_𝑗=1∫^𝑡𝑓₀Δ𝑇_𝑡𝑥_𝑗,𝑡².𝑑𝑡 (13) 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: 𝑈^𝑙+1−𝑈^𝑙𝑈^𝑙≤𝑒, (14) 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: ∇×𝐸=𝑖𝜔𝜇𝐻,∇⋅𝐸=0,∇⋅𝐻=0,∇×𝐻=−𝑖𝜔𝜀^∗𝐸, (15) 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 𝜀^∗=𝜀+𝑖𝜎𝜔. (16) The volumetric EM power deposition, which averaged over time, is calculated by 𝑄_𝑟=𝜎|𝐸|²2, (17) 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.

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 𝑌(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.

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

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 𝑈=200 W/m^2°C as well as the desired one, which is obtained from the application of a known 𝑈 at 𝑡=1000s. 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 𝑡=1000 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 𝑡=1000 s for two different thermal conductivities of casing layer. Solid line: 𝑘cl=0.14 W/m°C, 𝑈=175 W/m2°C; dashed line: 𝑘cl=0.11 W/m°C, 𝑈=150 W/m2°C.

Figure 8 

Temperature distribution at 𝑡=1000 s for two different heat convection coefficients of water. Solid line: ℎ=200 W/m2°C,𝑈=175 W/m2°C; dashed line: ℎ=350 W/m2°C, 𝑈=200 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, 𝑈_𝜀=𝑈+Δ𝑈. (A.1) Then, the temperature undergoes the variation Δ𝑇_𝑡(𝑥,𝑡) as follows: 𝑇_𝑡𝜀(𝑥,𝑡;𝑈)=𝑇_𝑡(𝑥,𝑡;𝑈)+Δ𝑇_𝑡(𝑥,𝑡;𝑈), (A.2) 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: 𝜌_𝑡𝑐_𝑡𝜕Δ𝑇_𝑡(𝑥,𝑡)𝜕𝑡=𝑘_𝑡𝜕²Δ𝑇_𝑡(𝑥,𝑡)𝜕𝑥²−𝜔_bl𝜌_bl𝑐_blΔ𝑇_𝑡(𝑥,𝑡)0<𝑥<𝑎,Δ𝑇_𝑡(𝑥,𝑡)=0𝑥=0,−𝑘_𝑡𝜕Δ𝑇_𝑡(𝑥,𝑡)𝑇𝜕𝑥=Δ𝑈_∞−𝑇_𝑡(𝑥,𝑡)−𝑈Δ𝑇_𝑡(𝑥,𝑡)𝑥=𝑎,Δ𝑇_𝑡(𝑥,𝑡)=00<𝑥<𝑎,𝑡=0. (A.3)(A.4)(A.5)(A.6) The second step is to calculate the variation of the objective function, Δ𝐺(𝑈), which is given by Δ𝐺(𝑈)=^𝑡𝑓₀^𝑎₀2Δ𝑇_𝑡𝑥_𝑗,𝑡𝑇_𝑡𝑥_𝑗,𝑡−𝑌_𝑡𝑥_𝑗,𝑡×𝛿𝑥−𝑥_𝑗𝑑𝑥𝑑𝑡. (A.7) We obtain the variation of the augmented functional, Δ𝐿(𝑇,Ψ_𝑡,𝑈), by considering the variation problem given by (A.3) and the variation of the objective function, (A.7), as follows: Δ𝐿=^𝑡𝑓₀^𝑎₀2Δ𝑇_𝑡𝑥_𝑗,𝑡𝑇_𝑡𝑥_𝑗,𝑡−𝑌_𝑡𝑥_𝑗,𝑡×𝛿𝑥−𝑥_𝑗+𝑑𝑥𝑑𝑡_{𝑡𝑓0𝑎0}Ψ_𝑡(𝑥,𝑡)𝑘_𝑡𝜕²Δ𝑇_𝑡(𝑥,𝑡)𝜕𝑥²−𝜌_𝑡𝑐_𝑡𝜕Δ𝑇_𝑡(𝑥,𝑡)𝜕𝑡−𝜔_bl𝜌_bl𝑐_blΔ𝑇_𝑡(𝑥,𝑡)𝑑𝑥𝑑𝑡. (A.8) Using integration by parts in (A.8), ^𝑡𝑓₀^𝑎₀𝜕²Δ𝑇_𝑡(𝑥,𝑡)𝜕𝑥²Ψ_𝑡()=𝑥,𝑡𝑑𝑥𝑑𝑡_𝑡𝑓0𝜕Δ𝑇_𝑡(𝑎,𝑡)Ψ𝜕𝑥^𝑡(𝑎,𝑡)−𝜕Ψ_𝑡(𝑎,𝑡)𝜕𝑥Δ𝑇_𝑡()−𝑎,𝑡𝜕Δ𝑇_𝑡(0,𝑡)Ψ𝜕𝑥^𝑡()+0,𝑡𝜕Ψ_𝑡(0,𝑡)𝜕𝑥Δ𝑇_𝑡(+0,𝑡)𝑑𝑡_{𝑡𝑓0𝑎0}𝜕²Ψ_𝑡(𝑥,𝑡)𝜕𝑥²Δ𝑇_𝑡(𝑥,𝑡)𝑑𝑥𝑑𝑡,_{𝑎0𝑡𝑓0}𝜕Δ𝑇_𝑡(𝑥,𝑡)Ψ𝜕𝑡^𝑡()=𝑥,𝑡𝑑𝑡𝑑𝑥_𝑎0Δ𝑇_𝑡𝑥,𝑡_𝑓Ψ_𝑡(𝑥,𝑡)−Δ𝑇_𝑡(𝑥,0)Ψ_𝑡(−𝑥,0)𝑑𝑥_{𝑎0𝑡𝑓0}Δ𝑇_𝑡(𝑥,𝑡)𝜕Ψ_𝑡(𝑥,𝑡)𝜕𝑥𝑑𝑡𝑑𝑥, (A.9)(A.10) and by using the boundary conditions of the direct problem in variations (Equations (A.3)–(A.6) and the rearrangement of the mathematical terms) the final form of (A.8) can be written as follows: Δ𝐿=^𝑡𝑓₀^𝑎₀𝑘_𝑡𝜕²Ψ_𝑡(𝑥,𝑡)𝜕𝑥²+𝜌_𝑡𝑐_𝑡𝜕Ψ_𝑡(𝑥,𝑡)𝜕𝑡−𝜔_bl𝜌_bl𝑐_blΨ_𝑡()𝑇𝑥,𝑡+2_𝑡𝑥_𝑗,𝑡−𝑌_𝑡𝑥_𝑗,𝑡×𝛿𝑥−𝑥_𝑗Δ𝑇_𝑡()−𝑥,𝑡𝑑𝑥𝑑𝑡_𝑡𝑓0𝑘_𝑡Ψ_𝑡(0,𝑡)𝜕Δ𝑇_𝑡(0,𝑡)+𝜕𝑥𝑑𝑡_𝑡𝑓0Ψ_𝑡(𝑎,𝑡)𝑈−𝑘_𝑡𝜕Ψ_𝑡(𝑎,𝑡)𝜕𝑥Δ𝑇_𝑡()−𝑎,𝑡𝑑𝑡_𝑎0𝜌_𝑡𝑐_𝑡Ψ_𝑡𝑥,𝑡_𝑓Δ𝑇_𝑡𝑥,𝑡_𝑓+𝑑𝑥_𝑡𝑓0Δ𝑈Ψ_𝑡(𝑎,𝑡)𝑇_𝑡(𝑎,𝑡)−𝑇_∞𝑑𝑡. (A.11)

To satisfy the optimality condition, that is, Δ𝐿(Δ𝑇,Ψ_𝑡,Δ𝑈)=0, it is required that all coefficients of Δ𝑇_𝑡(𝑥,𝑡) are to be vanished. Finally, the following term is left: Δ𝐿𝑇,Ψ_𝑡,𝑈=^𝑡𝑓₀Δ𝑈Ψ_𝑡(𝑎,𝑡)𝑇_𝑡(𝑎,𝑡)−𝑇_∞𝑑𝑡. (A.12) Since we know that 𝑇_𝑡(𝑥,𝑡) and Ψ_𝑡(𝑥,𝑡) are, respectively, the direct and the adjoint problem solutions, it follows that 𝐿𝑇,Ψ_𝑡,𝑈=𝐺(𝑈), (A.13) and then Δ𝐿𝑇,Ψ_𝑡,𝑈=Δ𝐺(𝑈). (A.14) Based on the definition, the variation of the functional 𝐺(𝑈) can be written as follows: Δ𝐺(𝑈)=⟨∇𝐺,Δ𝑈⟩, (A.15) where the associated inner scalar product ⟨,⟩ of two functions 𝑓(𝑡) and 𝑔(𝑡) in the Hilbert space of 𝐿₂(0,𝑡_𝑓) is defined by 𝑓(𝑡),𝑔(𝑡)_𝐿2=^𝑑_𝑐𝑓(𝑡)𝑔(𝑡)𝑑𝑡. (A.16) From the definition of the scalar product, the variation of the functional 𝐺(𝑈) can be written as follows: Δ𝐺(𝑈)=⟨∇𝐺(𝑈),Δ𝑈⟩_𝐿2=∇𝐺(𝑈)Δ𝑈. (A.17) A comparison between (A.12), (A.14), and (A.17) 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
bl:	Blood
𝑐:	Core
cl:	Casing layer
𝑓:	Final
𝑖:	Initial time
𝑡:	Tissue
∞:	Ambient.

Superscripts

𝑙:	Iteration number.