Online feedback focusing algorithm for hyperthermia cancer treatment

Abstract

Purpose: Magnetic resonance (MR) imaging is increasingly being utilized to visualize the 3D temperature distribution in patients during treatment with hyperthermia or thermal ablation therapy. The goal of this work is to lay the foundation for improving the localization of heat in tumors with an online focusing algorithm that uses MR images as feedback to iteratively steer and focus heat into the target.

Methods: The algorithm iteratively updates the model that quantifies the relationship between the source (antenna) settings and resulting tissue temperature distribution. At each step in the iterative process, optimal settings of power and relative phase of each antenna are computed to maximize averaged tumor temperature in the model. The MR-measured thermal distribution is then used to update/correct the model. This iterative procedure is repeated until convergence, i.e. until the model prediction and MR thermal image are in agreement. A human thigh tumor model heated in a 140 MHz four-antenna cylindrical mini-annular phased array is used for numerical validation of the proposed algorithm. Numerically simulated temperatures are used during the iterative process as surrogates for MR thermal images. Gaussian white noise with a standard deviation of 0.3°C and zero mean is added to simulate MRI measurement uncertainty. The algorithm is validated for cases where the source settings for the first iteration are based on erroneous models: (1) tissue property variability, (2) patient position mismatch, (3) a simple idealized patient model built from CT-based actual geometry, and (4) antenna excitation uncertainty due to load dependent impedance mismatch and antenna cross-coupling. Choices of starting heating vector are also validated.

Results: The algorithm successfully steers and focuses a tumor when there is no antenna excitation uncertainty. Temperature is raised to ≥43°C for more than about 90% of tumor volume, accompanied by less than about 20% of normal tissue volume being raised to a temperature ≥41°C. However, when there is antenna excitation uncertainty, about 40% to 80% of normal tissue volume is raised to a temperature ≥41°C. No significant tumor heating improvement is observed in all simulations after about 25 iteration steps.

Conclusions: A feedback control algorithm is presented and shown to be successful in iteratively improving the focus of tissue heating within a four-antenna cylindrical phased array hyperthermia applicator. This algorithm appears to be robust in the presence of errors in assumed tissue properties, including realistic deviations of tissue properties and patient position in applicator. Only moderate robustness was achieved in the presence of misaligned applicator/tumor positioning and antenna excitation errors resulting from load mismatch or antenna cross coupling.

• Phased array optimization
• automatic feedback control
• feedback focusing
• hyperthermia
• power deposition
• effective perfusion
• eigenvector

Introduction

Focused adjuvant hyperthermia that maintains a high temperature in tumor has been found to enhance the efficacy of radiotherapy and chemotherapy [1–5]. For example, hypoxic tumor regions may be re-oxygenated and thereby become radio-sensitiv [6], [7] when the tumor temperature is raised to mild levels (e.g. 40–44°C) for 30–60 minutes. Focused heating may also be used to trigger the local release of anti-cancer drugs encapsulated in thermally sensitive liposomes [8–10]. Such focused adjuvant hyperthermia facilitates selective drug release in tumor, thereby reducing normal organ toxicity [11]. There are also increasing reports that support the efficacy of stand alone heating at even higher temperatures (i.e. >48°C), such as thermal ablation therapy [12], [13].

Continued effort over the past two decades has produced devices that can deliver heat energy more controllably to the desired tumor target than previous equipment which had poor adjustability [14]. Despite the great potential of hyperthermia, delivering a safe and effective hyperthermia treatment remains a difficult task, due to the lack of a theoretically rigorous and experimentally validated bio-heat transfer model, complicated patient anatomical structures and geometries, and uncertainties associated with tissue properties [15] and blood perfusion, which vary significantly both spatially and temporally during treatm [16], [17]. However, adequate temperature elevation in the tumor is important for the success of hyperthermia. The presence of cold regions inside the tumor target has been linked to the failure of hyperthermia treatments [18–20]. Therefore, an online treatment-focusing algorithm, capable of dynamically adapting the model to uncertainties and changes in model parameters during treatment, would greatly enhance the clinical efficacy of hyperthermia.

Using MR to control hyperthermia treatments is attractive, especially in light of recent research reporting strong correlation between MR measured temperatures and hyperthermia treatment planning/outcomes [21–23]. As opposed to control based on a small number of measurement points from implanted temperature probes, an MR-based online controller can use information from the entire 3D temperature distribution as obtained from MR thermal images in multiple slices, to adjust the parameters of the heating device. Since more detailed spatial temperature information is used, better control is expected to overcome the temperature heterogeneity mentioned in the previous paragraph. MRI-based controllers have been reported in literature [24]. However, due to the complexity of the problem both temporally and spatially, few controllers operate in 3D [25], [26].

Controlling temperature in three dimensions implicitly relies on optimization of the temperature distribution, a subject that has been extensively studied in hyperthermia literature [27–43]. Boag and Leviatan [27] proposed an important analytical approach that reduces the computational burden of determining an optimal heating vector that simultaneously heats tumor and protects normal tissues by optimizing the ratio of tumor to normal tissue SAR. Later Kremer and Louis [28] gave a rigorous mathematical analysis of this analytical approach and Bohm et al. [29] presented a review of analytical perspectives of this approach. An important methodology extension was made by Das and Clegg et al. [30] to analytically optimize the ratio of tumor to normal tissues temperature so that blood perfusion effects are considered in hyperthermia optimization/control. Further advances in theoretical optimization have dealt with treatment practices [31], antenna array desi [32], [33], better design of objective functions or optimization methods that are more capable of associating practical optimal treatment to analytically/numerically optimized results [34–36], and reduction of computation loads using numerical techniques, e.g. combing quasi-static zooming [44] and optimization, to facilitate real-time controller design [37]. More recently, there have been reports on further extending the existing theoretical framework to transient temperature optimization/control [38–43].

The proposed work aims to overcome a limitation associated with the previous theoretical framework, i.e. the need for pretreatment model identification based on M² linearly independent predetermined heating vectors for an M-antenna phased array applicator. Power deposition during this process usually heats the patient inefficiently, since it is not designed to be focused at the tumor. Moreover, this process could be prohibitively long when there are many antennas. Our basic idea relies on quickly building a simplified patient model to provide a reliable first heating vector, and then iteratively updating our model in a least-square error (LSE) sense. Thus, each successive heating vector (at each iteration step) is incrementally closer to the true optimal heating vector, while potentially providing useful tumor heating from the very first iteration. In comparison, previous studies by Kowalski et al. [39], [40] relied on M² pretreatment model identification steps; more recently, they added a single-input-single-output proportional-integrate temporal power controller to enhance heating efficacy [40].

We propose a feedback-focusing algorithm for use prior to heat treatment, for initial planning optimization of antenna power and phase drive parameters. Thermal distributions obtained from MR images are fed back to update the model at every step in the iterative process until convergence (i.e. no significant updates in the driving vector). In this work, numerically simulated temperature distributions with some added Gaussian random error are used to mimic MR thermal images that would be obtained during an MR monitored hyperthermia treatment. The feedback-focusing algorithm is numerically simulated and assessed for the case of a human upper leg heated with a 140 MHz four-antenna cylindrical phased array hyperthermia applicator. The first driving vector (set of four phase and magnitude drive parameters) in the iterative feedback sequence is assumed to be obtained from an erroneous initial model. The initial model used for the simulations (prior to feedback-focusing) is assumed to be an incorrect characterization of the actual tissue target in four tests: (1) 30% deviation in the initial model electric properties and 3% deviation in the initial model thermal properties from that of the actual tissue target, (2) patient and tumor position in the initial model erroneously assumed to be shifted by 2 cm from the actual position, (3) a simplified homogeneous muscle equivalent tissue model with the same exterior geometry as the patient, and (4) patient position in the initial simple model erroneously assumed to be shifted by 2 cm from the actual position. Convergence of the algorithm to provide efficient heating of the target tumor while minimizing heating in surrounding non-tumor tissue is evaluated for each of these challenging cases.

Methods

Derivation of parametric model combining power and temperature

The mathematical models involved in simulating a hyperthermia treatment are the power model, which computes the power distribution in the heated domain, and the bio-heat transfer model, which uses the power distribution as input to compute the temperature distribution. This work assumes that the heating is caused by multiple electromagnetic (EM) sources. The EM wave model computes the electric fields which give rise to the power deposition field. The absorbed power deposition is an interior forcing source term in the bio-heat transfer equation. The EM field at any spatial location can be written as follows [45].where M is the number of antennas, E_i is the electric field from antenna i, and u_i is the complex source input value for antenna i. The above equation states that the resultant electric field is the weighted sum of the individual electric fields generated by each of the EM wave sources. The tissue-absorbed power distribution of this resultant electric field can then be written in Hermitian form [27].where superscript H denotes complex conjugate transpose. The exact form of the bio-heat transfer equation, the temperature distribution with zero boundary conditions can be expressed in terms of Green's function (the impulse response function, a result of applying both the superposition property and the homogeneity property) as follows [46].Here i, j and k are node indices along the x, y and z direction. The point-wise system matrices ([S_{i, j, k}]) are Hermitian, and is the driving vector. The temperature solution expressed by Equation 3 is valid for any linear model such as the k-effective model [47], the Pennes bio-heat transfer model [48], or a combination of both. For a heating system with M antennas/sources, the elements of the point-wise matrices [S_{i, j, k}] can be computed from either M² linearly independent MR thermal images or M² linearly independent bio-heat transfer simulations [30]. Following the determination of the point-wise matrices, they may be used to determine the temperature distribution for any arbitrary antenna settings.

The optimal driving vector that focuses power in the tumor is determined as the eigenvector of the total tumor system matrix corresponding to the maximum eigenvalue. The total tumor system matrix is given as follows.

Derivation of approximated combined model

The drawback of solving for M² unknowns in Equation 4 during the course of an actual treatment is that it requires thermal images corresponding to the application of M² linearly independent heat settings (starting from a uniform 37°C temperature distribution) [41–43]. This is not practical from the point of view of time, since it requires the thermal effect of each new heat setting to almost completely dissipate before imaging the effect of the next heat setting in the sequence. This imaging sequence is illustrated schematically in Figure 1(a). It is more practical to allow only a short cooling time between heat settings during the initial feedback power setting characterization studies, as illustrated in Figure 1(b). To mimic Figure 1(a), thermal images acquired at the more closely spaced time points shown in Figure 1(b) must be corrected by the feedback control algorithm to remove the cumulative effect of previous heat settings. The methodology used to remove the residual heat accumulated from previous settings is to fit an exponential fall-off to the cumulative temperature decay from the previous settings, an idea based on ‘effective perfusion’ (previously investigated and applied in other hyperthermia studies [49–55], [17], [56]. This is illustrated by the dashed line in Figure 1(b), wherein the temperature rise from the current setting is subtracted from the fitted decay of the previous heat settings (line X-Y in Figure 1(b)). The temperature image corresponding to this subtraction (for each pixel of the image) is used as the corrected temperature rise image for that setting. In the bio-heat transfer Equation [48]:where T, Q, ρ, C_t, C_b, k, ω_b, are the temperature, power deposition, tissue density, specific heat of tissue, specific heat of blood, thermal conduction, and perfusion, respectively. Exponential decay follows from the assumptions that (a) the perfusion term (ω_b) is constant with temperature, (b) the perfusion effect is much larger than the thermal conduction term (k), and (c) there is no input power (Q) [17], [49–55], [56]. The corrected thermal images can then be used as inputs to solve for the terms in [S_{i, j, k}] (Equation 4).

Figure 1. Sketch showing the superposition of a series of heat-up and cool-down processes for a particular point in the heated target. (a) Ideally, heat from the previous settings should almost completely dissipate before the next setting is applied. Images acquired at the time points shown by circles can be subtracted to directly obtain the temperature for each heat setting (denoted by line X-Y for the second heat setting). (b) In reality, due to time limitations, the effect of the previous settings is exponentially decayed. The image corresponding to the second heat setting is denoted by line X-Y.

Derivations of least-squares errors update algorithm

Completely solving for the terms in [S_{i, j, k}] requires M² linearly independent thermal images (obtained as described in the previous paragraph), which requires some initial time expenditure. To minimize this period of unfocused heating required for system characterization, rather than waiting to gather all M² thermal images the currently available thermal images (during the process of feedback) could be used to determine a least-squares (LS) error approximation to the desired system matrices [S_{i, j, k}]. More precisely, we start with an initial guess driving vector (driving vector that attempts to focus heat in the tumor) and obtain a corresponding thermal image, from which the first LS approximated system is determined. Using this LS approximated system matrix, the second driving vector is determined as the eigenvector corresponding to the maximum eigenvalue of the system matrix in Equation 4 (optimal driving vector). This process is repeated until all M² thermal images are gathered, at which point the system matrix is exactly determined, in theory. Each successive driving vector attempts to be useful by maximizing tumor heating. However, if the optimal driving vectors from one step to the next are very similar, the resulting thermal images might not be sufficiently linearly independent to robustly solve for [S_{i, j, k}]. Thus, greater than M² thermal images might be required to fully define the system matrix.

Numerical simulation setup

To validate the proposed feedback-focusing algorithm, numerical simulations were run for EM hyperthermia heating of a human upper leg (see Figure 2). The heating device is a simulated cylindrical array of four equispaced paired dipole antennas, which is used routinely in the Duke Radiation Oncology Clinic. The applicator is 23 cm in diameter, 27.5 cm in length, and driven at 140 MHz. It consists of eight copper foil strip dipole antennas connected in parallel pairs and printed on the inner surface of a cylindrical plastic shell. The radio frequency input to each antenna is across the center of the U-shaped paired dipole arms. For effective coupling to the patient, an expandable thin plastic membrane filled with water (water bolus) fills the space between the patient (e.g. leg and arm) and the cylinder. The human thigh model was approximately centered within this heating device. The driving vectors determined during the feedback update process were applied for two minutes to heat up the leg, followed by a one minute cool-down (see Figure 1(b)). These times for heat-up and cool-down were determined based on clinical experience and our previous study [57]. Simulated thermal images were used as surrogates for those that will be acquired in practice using, for example, MR thermal imaging. These images were computed in the heated volume at time points indicated by the circles in Figure 1(b) (representation of temperature-time history for a spatial point in the thermal image), corresponding to the end of heat-up and end of cool-down for each driving vector in the feedback update sequence. The steps involved in the computation that provided surrogate thermal images during the feedback update process were as follows: (1) Finite Element (FE) software (HFSS -Ansoft Corporation (Pittsburgh, PA)) was used to compute the electric field distribution in the leg and output the fields in a uniform spatial grid of 4 mm, (2) the electric fields were used as inputs to an in-house developed finite difference software to compute the thermal distribution. A constant uniform water bolus temperature of 37°C was used as the boundary condition. Numerical grid sizes of 0.5 s and 4 mm were used for temporal and spatial discretization, respectively. The property values used are summarized in Table I.

Figure 2. Patient upper leg finite element model, shown within the 140 MHz mini-annular phased array heating apparatus (red for tumor, blue for bone, green for muscle and yellow for fat).

Table I.  Nominal properties’ values (assumed for 140 MHz).

	Relative electric permittivity, εR, εR = ε/ε0	Electrical conductivity,  σ (mho/m)	Density,  ρ (kg/m^3)	Specific heat, Ct (J/kg/°K)	Thermal conductivity k (W/m/°C)	Blood perfusion wb (kg/m^3/s)
Blood	–	–	–	3770	–	–
Bone	10	0.05	1600	1300	0.40	0.0
Fat	10	0.06	880	2387	0.23	1.1
Muscle	75	0.75	1050	3639	0.50	3.6
Tumor	59	0.65	1050	3639	0.50	1.8
Water	81	0.0002	996	4178	0.61	–

Permittivity of vacuum, ε₀ = 8.854 * 10⁻¹² Farad/m and permeability of vacuum, μ₀ = 4 * π * 10⁻⁷ Henry/m. The permeabilities of all tissues were assumed to be equal to that of vacuum.

Robustness studies

The first optimal driving vector in the feedback-focusing sequence is based on an initial model of the heated target. The greater this model is in error from the actual situation, the more likely it is that the feedback control procedure will fail to converge. The effects of various error sources are investigated to validate the robustness of the proposed algorithm in finding an appropriate solution.

Tissue property variability

The first source of error considered was incorrect estimation of tissue properties. Although experimental measurements of electric and thermal properties of human tissues are documented in the literature, actual values vary considerably between patients [15], across tissue types [58], and as a function of time during treatment. To account for these uncertainties, initial values of patient tissue properties such as electrical conductivity and permittivity, were assumed to be deviated ±30% from the actual properties. Also, thermal conductivity, density and blood perfusion were assumed to be deviated ±3% from the actual properties. Patient property deviations in the initial model are summarized in Table II.

Table II.  Property deviations.

	Relative electric permittivity, εR, εR = ε/ε0	Electrical conductivity,  σ (mho/m)	Density,  ρ (kg/m^3)	Specific heat, Ct (J/kg/°K)	Thermal conductivity k (W/m/°C)	Blood perfusion wb (kg/m^3/s)
Blood	–	–	–	+3%	–	–
Bone	+30%	+30%	−3%	+3%	+3%	+3%
Fat	−30%	−30%	+3%	−3%	+3%	+3%
Muscle	+30%	+30%	−3%	+3%	+3%	−3%
Tumor	−30%	−30%	+3%	−3%	+3%	+3%

Patient position mismatch

The second source of error considered was patient positioning. Even though techniques of patient positioning and tumor targeting have improved with technology (e.g. with better MRI guidance), it is difficult to position the patient precisely within the flexible bolus. Thus, unlike radiation therapy, where better than 0.5 cm in patient positioning is expected, more significant patient position error is unavoidable in a hyperthermia treatment [32]. To characterize the effect of such an error, the actual patient position was assumed to be shifted by 2 cm in the +y direction (see Figure 2) from the position assumed in the initial model. Moreover, since building a detailed patient model can require 1 to 2 days if built manually [59] or ∼4 h if built semi-automatically [59], [60], a simple idealized model was considered instead. The simple model consisted of a homogeneous, muscle equivalent tissue volume, with the same external geometry as the actual patient leg.

Noise in measured temperature images

The third source of error considered was noise in the thermal images. Since the proposed algorithm is expected to use magnetic resonance thermal images (MRTI) as feedback to steer and focus power, a certain amount of noise is likely to be present. Therefore, to determine robustness of the controller in the presence of noise, Gaussian white noise with a standard deviation of 0.3°C was added to the simulated temperature field. This MRI standard deviation was determined from our recent phantom experiment, using the temperature-phase conversion coefficient determined in a previous study [61].

Antenna excitation uncertainty

The fourth source of error considered was antenna-load mismatch that produces discrepancies between the applied and actually obtained antenna magnitudes and phases. These discrepancies could be a result of equipment calibration changes over time or, more likely, load-dependent impedance mismatch and antenna cross-coupling. For example, the mini-annular phased array extremity heating applicator being used at our institution has been shown to have load-dependent deviations up to 15–20% in antenna phases and magnitudes. Table III summarizes the four sets of antenna discrepancies investigated here. For any particular set (table row), the actual antenna magnitudes were assumed to be scaled up/down from the desired magnitudes by the scale factors shown in the second column; the actual phases were assumed to be shifted from the desired values by the values shown in the third column. Since the tumor location is not at the center of the leg model (see Figure 1), the first set of discrepancies tends to increase tumor heating (since Ant. 1 and Ant. 2, closer to the tumor, have scaled up magnitudes), while the third set tends to de-focus tumor heating (since Ant. 3 and Ant. 4, away from the tumor, have scaled up magnitudes). The other two sets fall between these extremes. The simulations for antenna discrepancy used the deviated patient model. Both minimum and maximum eigenvector of the initial model were used as the first heating vector to test convergence.

Table III.  Uncertainty in antenna phase and magnitude settings.

	Magnitude discrepancy (scaling)				Phase discrepancy (degree)
Error condition	Ant. 1	Ant. 2	Ant. 3	Ant. 4	Ant. 1	Ant. 2	Ant. 3	Ant. 4
Set 1	1.2	2.0	0.5	0.8	20	30	−50	−10
Set 2	2.0	0.5	0.8	1.2	30	−50	−10	20
Set 3	0.5	0.8	1.2	2.0	−50	−10	20	30
Set 4	0.8	1.2	2.0	0.5	−10	20	30	−50
	Bottom	Left	Top	Right	Bottom	Left	Top	Right

Choices of starting heating vector

Lastly, to gauge robustness of the feedback-focusing algorithm to the choice of starting point, the first driving vector was set as the source combination that yielded the best or worst tumor heating in the initial model. Essentially, these settings were proportional to the eigenvector corresponding to the largest and smallest eigenvalue. Simulations were performed for the four different initial models: deviated standard model, shifted standard model, idealized model, and shifted idealized model. Each simulation involves 48 updates of the model using the simulated MR Images (2-min heat-up followed by a 1-min cool-down). Even though only 16 iteration steps are required, in theory, to achieve focus in the absence of errors, we extend iteration steps to observe if further improvement in focusing can be made and to see how long it would take to reach a steady state treatment and to heat enough tumor volume. This information can be used for further controller design if needed. MRI noise was included in all simulations.

Results

Simulations with no antenna discrepancy

To track heat steering and focusing in the tumor, the temperature distribution in the central slice (z = 0 cm) is shown in Figures 3 and 4 for an increasing number of feedback iterations. Figure 3 illustrates the iteration history when the first heating vector was determined as the maximum eigenvector of Equation 4 (i.e. best initial guess of four antenna phase and magnitude values) of the initial model with a ±30% electric property deviation and ±3% thermal property deviation (see Table II for details). The temperature distribution in the tumor initially appeared to be focused (first two steps), but subsequently lost focus starting from the fourth to the tenth step, possibly as a result of noise and approximations used to reconstruct the model (see Methods). However, as the number of iterations increased to close to the full reconstruction at the 15th to 16th step (4² = M² iterations), the feedback-focusing algorithm steered the focus back to tumor. Though a slight de-focus appeared again from the 16th to 20th step, this time power deposition was confined within a smaller region closer to tumor and was finally steered back to tumor position by the 21st step. Subsequently, the focus remained centered in tumor with little change (steps 24–48).

Figure 3. Temperature distribution at a central axial slice through tumor (z = 0 cm) at the end of the feedback iteration number indicated in top right of each panel. The first driving vector was the maximum eigenvector from an initial model with property deviations. The outermost ring is fat, within which is the muscle region. The larger irregular object on the left side of the muscle region is tumor, and the smaller irregular object at the center is bone. The colored regions are: red for ≥43°C, yellow for ≥ 41°C, light blue for ≥40°C and blue for ≥39°C.

Figure 4. Temperature distribution at a central axial slice through tumor (z = 0 cm) at the end of the feedback iteration number indicated in top right of each panel. The first driving vector was the maximum eigenvector of the idealized model shifted by 2 cm. Geometry and labels are the same as Figure 3.

Similarly, Figure 4 illustrates the iteration history when the first heating vector was the maximum eigenvector of the initial idealized model with a 2 cm shift in the y direction (see Methods). Again, the initial focus was good but was sustained only for the first two iterations before losing focus for the fourth to twelfth steps. The algorithm steered the focus back to the tumor starting from the 13th step. Since the initial model involved greater error (simplified and purposely shifted), the algorithm required more iteration steps (23) to fully re-focus and reach a stable focus. However, a larger region of normal tissue at the top of leg, which was incidentally heated during the period of lost focus in the fourth to twelfth steps, was maintained at a higher temperature than desired (as seen in the anterior portion of the cross-section in Figure 4). This indicates a broader focus, but one that was still centered in the tumor target.

To visualize focusing in three dimensions, Figure 5 summarizes the temperature distributions corresponding to iteration 49 at three axial positions along the leg (z = 0 and ±3.6 cm), for the four different initial models. Figures 5(a) and 5(b) are for the cases where the first driving vector corresponds to the maximum eigenvector (good initial starting point) and minimum eigenvector (bad initial starting point), respectively. The results indicate that focusing was successful in all three axial locations along the leg when starting with the maximum eigenvector as driving vector (Figure 5(a)). Also, less normal tissue was detrimentally heated (i.e. temperature ≥41°C) when there were no erroneous shifts of leg position from the pre-planned location within the cylindrical array. In contrast, when a bad initial guess (minimum eigenvector) was used for the treatment planning, the algorithm was not able to correct for the shifted patient position or the imperfect ‘simplified’ geometry, as evidenced by the poor localization of heating in tumor seen in Figure 5(b). Figure 6(a) and (b) show the percentages of heated tumor volume and normal tissues volume for the four different initial models, for initial driving vectors corresponding to the maximum and minimum eigenvectors, respectively. Here, a tumor point is assumed to be therapeutically heated if its temperature is ≥43°C while normal tissue is defined as being detrimentally heated if ≥41°C. Note that, as a result of the physics underlying the heat dissipation process, normal tissues in the immediate vicinity of heated tumor will necessarily be at a higher temperature. Thus, those normal tissues falling within 1 cm around the tumor were not considered when determining the percentage of normal tissue heated to ≥41°C. For the case shown in Figure 6(a) of initial driving vector corresponding to the maximum eigenvector (i.e. good starting point), more than about 90% of tumor temperature was raised to ≥43°C for all four simulated initial models. Approximately 12–20% of normal tissue was raised to ≥41°C (except for the simplified tissue model with 2 cm shift in the y direction, where it was ∼25%). For the case shown in Figure 6(a) with initial driving vector derived from the minimum eigenvector (bad starting point), heat localization to the tumor was much worse when there was a 2 cm shift in the y direction. The percentage of tumor with temperature ≥43°C for these two initial models (shifted, or simplified and shifted model) could not be raised above 60%, with a corresponding expense of ∼20% of normal tissue temperature ≥41°C. If no position shift was included in the initial model, ∼90% of tumor was raised to ≥43°C. However, for the idealized initial model, ∼35% of normal tissue temperature was also raised to ≥41°C.

Figure 5. (a) Temperature distribution at different axial locations along the thigh at the 49th iteration, for an initial driving vector corresponding to the maximum eigenvector of the erroneous models described at the top of each panel. Geometry and labels are the same as Figure 3. (b) Temperature distributions at different z-locations using the 49th driving vector when the first heating vector was chosen as the minimum eigenvector for the four numerical validations. Geometry and labels are the same as in Figure 3. For patient position shift, unlike Figure 5(a), focusing is not centered within the tumor target.

Figure 6. (a) The percentage of heated tumor volume (T ≥ 43°C) and normal tissues (T ≥ 41°C), when initial driving vector is the maximum eigenvector of the five initial models. (b) For patient position shift, using the minimum eigenvector as the first heating vector only raises ∼60% tumor temperature above 43°C. (b) The percentage of heated tumor volume (T ≥ 43°C) and normal tissues (T ≥ 41°C), when initial driving vector is the minimum eigenvector of the five initial models. For patient position shift, using the minimum eigenvector as the first heating vector only raises ∼60% tumor temperature above 43°C.

Figure 7 measures the similarity of the applied heating vector at each iteration step to the best possible heating vector obtained for the standard model (no property deviation, applicator position shift, or model simplification). The absolute value of the dot product between the applied heating vector and the maximum eigenvector of the standard model is used as a measure of similarity (0 = no similarity, 1 = exactly similar). Figure 7(a) and 7(b) correspond to initial driving vectors proportional to the maximum and minimum eigenvectors. In Figure 7(a), at approximately the 25th iteration, the feedback algorithm steered to a stable vector that is close to the true optimal focusing vector, the maximum eigenvector of the standard model. A similar iteration history is given in Figure 7(b) when the first heating vector was determined as the minimum eigenvector for the four numerical validations. The absolute value of the dot product after convergence is summarized in Table IV.

Figure 7. (a) The dot product between the heating vectors at each feedback iteration step to the maximum eigenvector of the standard model, when the first heating vector was obtained from maximal eigenvector for the four numerical validations. (b) The dot product between the heating vectors at each feedback iteration step to the maximum eigenvector of the standard model, when the first heating vector was obtained from the minimum eigenvector for the four numerical validations.

Table IV.  List of absolute value of the dot product after convergence.

	Deviated standard model	Shifted standard model	Idealized model	Shifted idealized model
Max e-vector for the first heating	0.79	0.83	0.78	0.90
Min e-vector for the first heating	0.84	0.55	0.77	0.58

Simulations with antenna discrepancy

Figure 8(a) and (b) show the percentages of tumor volume heated to ≥43°C and normal tissue volume ≥41°C, when the initial driving vector is the minimum/maximum eigenvector of the deviated patient model, for the four sets of antenna excitation errors (Table III). Results in Figure 8(a) show that the first and fourth sets of discrepancies provide a similar focused heating capacity in terms of percentages of therapeutically heated tumor volume and detrimentally heated normal tissue volume. In contrast, the second and third sets are similar to each other and take more iteration steps to produce effective localization in tumor. Though all simulated situations have more than about 80% of tumor ≥43°C by the 49th iteration, the percentages of normal tissue volume ≥41°C are all more than about 30%.

Figure 8. (a) The percentage of heated tumor volume (T ≥ 43°C) and normal tissues (T ≥ 41°C), when initial driving vector is the maximum eigenvector from an initial model with property deviations in combination with four sets of antenna discrepancies (Table III). In general, more normal tissues are detrimentally heated in the presence of antenna discrepancy. (b) The percentage of heated tumor volume (T ≥ 43°C) and normal tissues (T ≥ 41°C), when initial driving vector is the minimum eigenvector from an initial model with property deviations in combination with four sets of antenna discrepancies (Table III). In general, more normal tissues are detrimentally heated in the presence of antenna discrepancy.

Discussion

The feedback-focusing algorithm presented here provides a solution to the important hyperthermia treatment problem of efficiently steering and focusing externally radiated radiofrequency energy into a tumor target at depth in the body. This algorithm is numerically validated to be robust to two common error sources: tissue property variability [15] and patient-applicator position mismatch relative to the treatment plan, as shown in Figures 3 and 5. Furthermore, the feedback-focusing algorithm appears to be robust even when the first driving vector is determined by maximizing the averaged tumor temperature of a simple simplified model with/without a 2 cm position mismatch. As shown in Figures 4 and 5A, the initial focus was lost in the first few iteration steps and then regained in subsequent iterations. The use of an idealized initial model relieves the burden of building an accurate and detailed patient model (1 to 2 days if built manually [59] and ∼4 h. if built semi-automatically [59], [60]) and the associated time-intensive numerical simulations. Depending on the initial model used, when the maximum eigenvector was used as the first heating vector Figure 6(a), the feedback-focusing algorithm successfully raised approximately 90–99% of the tumor temperature ≥ 43°C and only about 25% of normal tissue temperature ≥ 41°C. Therefore, the proposed feedback-focusing algorithm relieves the burden of building an accurate patient model and the associated time-intensive numerical computation. It is also robust to patient positioning errors (+/−2 cm) and variable patient electrical/thermal properties.

Feedback-focusing simulations were also conducted when the first driving vector was chosen to be the eigenvector corresponding to the smallest eigenvalue (bad initial guess) of the initial model. As shown in Figure 5(b), in the absence of 2 cm position mismatch in the initial model (standard or simplified), the feedback-focusing algorithm was still successful in steering and focusing power into the target region. However, when the initial model involved position mismatch, using the minimum eigenvector does not appear to provide enough information to allow the proposed algorithm to re-focus heating at the target region. As a result (see Figure 6(b)), only ∼60% of tumor volume was heated to a temperature ≥43°C, while ∼20% of normal tissue was heated to temperatures ≥41°C.

Lastly, this feedback-focusing algorithm demonstrates only a moderate robustness when there are errors in the source outputs. As shown in Figure 8, except for the third set of discrepancies, the algorithm still steered and focused power to heat at least ∼80% tumor volume therapeutically. However, at least ∼30% of the normal tissue volume was overheated. Therefore, the present algorithm requires some improvement to cope with antenna excitation errors that can occur as a result of load mismatch and antenna cross-coupling.

Our current work aims to provide useful tumor heating from the very start. Therefore, in contrast to other work [39], [40], system identification time required for model convergence is potentially well utilized. In Kowalski's algorithm [39], ∼50% of tumor volume was heated to a temperature ≥43°C, while ∼40% of normal tissue was heated to temperatures ≥41°C (using a higher frequency applicator with more antennas than that used in the current work). By comparison, our algorithm provided almost complete tumor heating with only ∼20% of normal tissues heated. This would appear to suggest that our algorithm has good potential for optimally focusing heat.

Further speed-up in model convergence time may be possible if one were able to obtain good quality thermal images in a shorter period of time, i.e. shorter time duration for each heat up iteration. MRTI techniques currently under development appear to be promising in this regard [25], [26], [62–64]. These techniques can potentially greatly reduce system identification time required to focus heating in the tumor, thereby not adding significantly to hyperthermia treatment time.

Five future research directions (the first three for EM synthesis and the remaining 2 for the bio-heat transfer process) can be incorporated to further test the effectiveness of this feedback process (i.e. increase the percentage of heated tumor volume and decrease the percentage of heated normal tissue): (1) use antennae operating at higher frequency [33], (2) use phased arrays with more rings of antennae along the z-directio[33], [65], [66] (e.g. the BSD-2000 Sigma-Eye system, which consists of three rings of dipoles along the longitudinal direction and four dipole antennae in each ring to allow longitudinal focusi [67], [68]), (3) use phased arrays with mo [33], [66] for better transversal focusing ability, (4) increase skin cooling [69], and (5) optimize the heat-up and cool-down sequence. According to the 2001 study made by Kroeze et al. [33], the first two methods provide more significant improvements in power focusing than the third method, since a finer spatial resolution is provided by using higher operating frequency, and since the tumor extends longitudinally and thus needs longitudinal focusing ability. However, heavier computational loads are involved, which makes this research topic more challenging. Regarding the two enhancements related to the bio-heat transfer process, a cooler boundary condition is definitely helpful, but cooling only affects tissues close to the skin surface [70]. Optimizing the heat-up and cool-down sequence implies greater accuracy of the MR thermal images, and thereby faster convergence of the feedback process. Further studies incorporating these enhancements are required to more fully understand their implications.

Conclusions

An iterative feedback-focusing algorithm for pretreatment planning optimization of heating from multi-antenna phased array hyperthermia applicators is shown to successfully steer and focus power into a target tumor for four typical challenging planning situations. Initial antenna drive vectors determined using either a detailed time-consuming patient specific tissue model or an easy-to-build simplified model, are shown to be robust starting points for the iterative treatment planning algorithm. Results demonstrate that the most fault tolerant focusing capability is possible when the initial phase and magnitude parameter set is chosen to be close to the correct set needed for appropriate tumor focus. But they also demonstrate an ability to return heating focus to tumor after about 25 iterations of the feedback algorithm starting with bad antenna parameter sets (derived from incorrect assumptions of tissue properties and tumor locations in the applicator). The proposed controller appears to be a good first step treatment planning approach to improve future MR-guided hyperthermia treatments with phased array applicators.

Acknowledgements

The authors would like to thank Dr Oana Craciunescu, for providing the patient CT images. The authors would also like to thank Ansoft Corporation (Pittsburgh, PA) for provision of software. The authors thank Dr James MacFall for comments on MRTI. The authors thank the reviewers for their suggestions and comments in improving this manuscript. The authors also thank Paul Turner from BSD for consultation on antenna calibration. This work was supported by grants NCI P01—CA042745-19 (SD) and NCI CA42745 (MWD).