Visualization of thermal washout due to spatiotemporally heterogenous perfusion in the application of a model-based control algorithm for MR-HIFU mediated hyperthermia

Abstract

Purpose

This article will report results from the in-vivo application of a previously published model-predictive control algorithm for MR-HIFU hyperthermia. The purpose of the investigation was to test the controller’s in-vivo performance and behavior in the presence of heterogeneous perfusion.

Materials and methods

Hyperthermia at 42°C was induced and maintained for up to 30 min in a circular section of a thermometry slice in the biceps femoris of German landrace pigs (n=5) using a commercial MR-HIFU system and a recently developed MPC algorithm. The heating power allocation was correlated with heat sink maps and contrast-enhanced MRI images. The temporal change in perfusion was estimated based on the power required to maintain hyperthermia.

Results

The controller performed well throughout the treatments with an absolute average tracking error of 0.27 ± 0.15 °C and an average difference of 1.25 ± 0.22 °C between $T_{10}$ and $T_{90}.$ The MPC algorithm allocates additional heating power to sub-volumes with elevated heat sink effects, which are colocalized with blood vessels visible on contrast-enhanced MRI. The perfusion appeared to have increased by at least a factor of ∼1.86 on average.

Conclusions

The MPC controller generates temperature distributions with a narrow spectrum of voxel temperatures inside the target ROI despite the presence of spatiotemporally heterogeneous perfusion due to the rapid thermometry feedback available with MR-HIFU and the flexible allocation of heating power. The visualization of spatiotemporally heterogeneous perfusion presents new research opportunities for the investigation of stimulated perfusion in hypoxic tumor regions.

Keywords:

• High intensity focused ultrasound
• Noninvasive thermometry
• Modeling
• Control systems engineering/treatment optimization
• Physiological effects of hyperthermia

Introduction

Mild hyperthermia, that is, the heating of malignant tissue to temperatures between 40 and 43 °C, sensitizes tumors for chemo- and radiotherapy (RT and CT) without introducing additional systemic toxicity [1,2]. This thermal enhancement effect has already been demonstrated in several clinical trials targeting a wide variety of tumors [3–9]. It arises from direct denaturation of proteins [10], the inhibition of repair mechanisms [11–13], and the killing of cells that are RT- and CT-resistant, but vulnerable to heat, that is, cells in the ‘S-phase’ [14,15] or in environments that are hypoxic and acidic due to inadequate perfusion [16–20]. The vasodilation following hyperthermia can also lead to increased perfusion in the tumor, resulting in sensitization to RT and CT via reoxygenation [21,22] and enhanced drug availability [23–26]. However, maintaining a high degree of control over the spatial temperature distribution during hyperthermia treatments is vital, as the clinical outcome would be negatively impacted by both the delivery of insufficient thermal dose [27–30] and vascular shutdown due to overheating [31].

Devices used for clinical hyperthermia applications are currently almost exclusively based on extracorporeal antenna systems emitting electromagnetic waves (EMW) in the radiofrequency and microwave range. These devices enable the application of energy across regions with diameters of tens of centimeters (regional hyperthermia, RHT), allowing to heat large tumor volumes. However, radiative phased array systems have limited capacity for the compensation of inhomogeneous heating rates and local heat sinks due to the lack of spatially resolved real-time feedback and the inherent lower limit to the size of the focus [32–35]. Temperature control is achieved via a-priori treatment planning with the help of dedicated model-based tools [33] and the software-assisted adjustment of treatment parameters [36] based on monitoring of the target tissue’s temperature via thermocouples [37] as well as patient feedback on off-target heating [38]. In the future, the use of magnetic resonance imaging (MRI)-based thermometry during EMW-hyperthermia may further improve clinical outcomes by providing the operator with 3D temperature information at higher spatial resolutions [39,40]. The systematic validation of a device for EMW-mediated RHT under MRI guidance is currently ongoing [34,35,41].

Magnetic Resonance-guided High-Intensity Focused Ultrasound (MR-HIFU) presents an alternative means by which hyperthermia can be applied. HIFU-mediated hyperthermia concepts initially also relied on a-priori planning of the treatment parameters [42,43]. However, after the development of fast proton resonance frequency shift (PRFS)-based thermometry [44,45], closed-loop temperature control concepts for MR-HIFU were developed quickly [46–49] and were applied in preclinical research on local drug delivery from temperature-sensitive liposomes [50–53]. Clinical HIFU transducers generate coherent ultrasound fields at single frequencies in the range of 0.5 ∼ 8 MHz [54]. The resulting focal heating volumes are only a few millimeters in diameter, as opposed to the volumes with diameters of tens of centimeters that are heated simultaneously with EMW devices. The heated volume in HIFU treatments can instead be increased by moving the focus, either via mechanical repositioning of the transducer (mechanical beam steering) or via coordinated phase shifts among the piezo elements on phased array transducers (electronic beam steering) [55–57]. This technique has already been applied clinically by Chu et al. and Heijman et al. for hyperthermia treatments of recurrent rectal cancer and dedifferentiated liposarcoma, respectively [58,59]. However, earlier calculations already showed that the high spatial control in power deposition of HIFU allows compensating for differences in local tissue cooling due to tissue perfusion and around thermally significant vessels, which is not feasible using a binary control algorithm together with electronic beams steering along fixed trajectories [60]. The logical next step in the development of HIFU technology for the application in hyperthermia is therefore to fully leverage the precision of HIFU and spatially resolved MR-thermometry by a transition from binary feedback control to a controller that calculates optimal sonication patterns. Model predictive control (MPC) is a modern control scheme that is widely used in the industry [61,62] and has already been applied in algorithms for thermal ablation therapy [63–65]. MPC algorithms use a model of the system’s dynamic behavior (the state-space model) to predict the state of the system at a future time point in response to the available control inputs. Using a cost function to define the target state of the system, the best control action at a given time point is found via the variation of the control inputs such that the cost function assumes its minimum and all state- and input variables satisfy their respective constraints. This allows the processing of multidimensional inputs like temperature maps into multidimensional outputs like heating power distributions, thereby enabling the online adaption of sonication patterns to the current distribution and strength of heat sinks.

While MR-HIFU mediated hyperthermia has not yet been widely adopted clinically alongside EMW-mediated RHT, the combination of rapid PRFS-thermometry and precise heating provides the potential for unmatched temporal and spatial temperature control. Recently we presented the architecture of a model predictive control (MPC) algorithm for MR-HIFU mediated local hyperthermia and characterized its performance in a perfused tissue-mimicking phantom in comparison to the current state-of-the-art binary feedback controller [57,66]. We showed that the MPC algorithm yields superior performance in-vitro via the targeted compensation of heat loss due to diffusion and local heat loss induced by an artificial blood vessel. In this study, we present results from the in-vivo application of the MPC algorithm in porcine thigh muscle. The aim is to evaluate the controller’s performance in live tissue and to investigate the interplay between its flexible heating power allocation and spatiotemporally varying perfusion.

Materials and methods

Animal model

The performance of an MPC algorithm for the creation and maintenance of local hyperthermia was investigated in muscle tissue of the biceps femoris of five healthy, female German landrace pigs between 35 and 39.5 kg in weight. The experimental protocol was approved by the State Agency for Nature, Environment and Consumer Protection of North-Rhine Westphalia, Germany.

Starting 12 h prior to the experiment, the pigs were fasted with free access to water. At the beginning of the experiment, the animals received an intramuscular injection of atropine (0.02 mg/kg; WDT, Germany), tiletamin/zolazepam (10 mg/kg; Zoletil, Virbac, Germany), and xylazine (2 mg/kg; Rompun 2%; Bayer; Germany). Anesthesia was induced and maintained intravenously using propofol (induction: 1.66 mg/kg i.v., maintenance: 4.0–9.5 mg/kg/h as required; Fresenius, Bad Homburg, Germany). Analgesia was achieved by intravenous administration of buprenorphine (0.02 mg/kg; Buprenovet Multidose, Bayer, Germany). The pigs were intubated and ventilated under pressure control (Hamilton C1, Heinen + Löwenstein, Switzerland) at 30% oxygen, 14 breaths/min, and positive end-expiratory pressure (PEEP) of 5–8 mmHg. The tidal volume was adjusted as needed to maintain normocapnia (PaCO₂ 35–45 mmHg). A catheter with three lumens (18 G; Arrow International, Reading, USA) for the continuous administration of propofol and Ringer’s solution (5 ml/kg/h; Fresenius Kabi, Germany) and the administration of other medication was placed in the right external jugular vein using the Seldinger technique.

The animals were placed inside the MRI scanner in right lateral recumbency with the biceps femoris above the acoustic window of the HIFU table (Figure 1(A)). A rubber tube containing MR-visible oil capsules was affixed to the right thigh to allow localization of the targeted volumes during dissection (Figure 1(B)). Acoustic contact was achieved by shaving, the application of degassed ultrasound gel and the interposition of an acoustic coupling gel pad (Aquaflex^® Ultrasound Gel Pad 4 × 27.5 × 27.5 cm, Parker Laboratories, Farfield, NJ, USA). Degassed and deionized water was used to ensure acoustic coupling between the acoustic membrane of the HIFU table and the gel pad. After positioning of the animal in the scanner was completed, a muscle relaxant (initial bolus: 0.1 mg/kg, infusion: 36 mg/h/kg; Pancuroniumbromid Rotexmedica 2 mg/ml, Panpharma GmbH, Trittau, Germany) was administered to prevent shivering.

Figure 1. Experimental setup. (A) The animal is positioned on the HIFU table such that the biceps femoris can be sonicated. (B) The rubber tube containing fish oil capsules is visible in MRI and is used to document the target site. (C) MRI with temperature map overlay. Yellow contour: Drift correction ROI with manually excluded areas. Red contour: Target ROI.

Up to three hyperthermia treatments with a target temperature of 42 °C were performed in each animal and were each continued for up to 30 min. In total, n = 10 treatments were performed throughout the study. Each treatment targeted a different point in the biceps femoris of the right leg. The target points were always separated by at least 3 cm. After the end of each treatment, temperature mapping was continued for up to 15 min and a minimum waiting time of 30 min was maintained between sonications to ensure the restoration of thermal equilibrium. The target points were positioned in the highly perfused areas close to the femur to maximize the likelihood of blood vessels being present in the vicinity (Figure 1(C)). The animals were returned to the stable eight hours after the beginning of anesthesia and were allowed to recover for 48 h under observation. Finally, the animals were euthanized using sodium-pentobarbital (150 mg/kg i.v., Euthadorm^®, CP-Pharma, Burgdorf, Germany).

Experimental setup

The hyperthermia treatments were performed on a clinical MR-HIFU system (3T Achieva^®, Philips Healthcare, Best, The Netherlands, and Sonalleve^® V2 HIFU, Profound Medical, Mississauga, Canada). During single-point sonications at 1.2 MHz in water, the transducer in this system generates ellipsoid focal volumes with 6 dB diameters of 1.7 mm and 13.6 mm perpendicular and along the transducer axis, respectively. During such sonications, 48.8% of the transducer’s output power pass through the focal volume and a time-averaged peak focal intensity $I_{sp}$ of 57.5 W/cm² per Watt of nominal output power is reached. Treatment planning and feedback control of the temperature during treatment was performed using a python-based software developed in-house and a model predictive control algorithm. The python-based software provides a user interface with tools for the planning and monitoring of treatments and communicates with the hardware using the pyMRI and pyHIFU packages of Zaprozan et al. [67].

Magnetic resonance imaging

All images were acquired using the HIFU table’s window coil and the HIFU pelvis coil (Model 905051-F, Philips Healthcare, Best, The Netherlands) on a 3T MRI scanner (3 T Achieva^®, Philips Healthcare, Best, The Netherlands).

Pretreatment imaging

The treatment planning was performed using a 3D T1-weighted fast field echo (FFE) without fat suppression. The used parameters were TE = 2.3 ms, TR = 6.1 ms, FA = 20°, acquisition matrix = 88 × 166, number of slices = 167, reconstructed (acquired) voxel dimensions = 0.49 × 0.49 × 1.50 mm³ (1.52 × 1.53 × 3.00 mm³), FOV = 140 × 261 × 250 mm³, fold-over direction = RL. Constant Level Appearance (CLEAR) was used to achieve higher signal homogeneity and Sensitivity Encoding (SENSE, P reduction = 1.2, S reduction = 1) was used to accelerate imaging. This resulted in an acquisition duration of 5.5 min.

MR thermometry

The imaging sequence used to monitor the treatment progress was an RF-spoiled gradient echo sequence with TR = 30 ms, TE= 19.5 ms, FA = 19.5°, FOV = 288 × 288 mm³, 6 slices, and voxel dimensions = 1.8 × 1.8 × 7 mm³, resulting in an acquisition time of 3.7 s for a full set of slices and a resolution of 160 × 160 voxels per slice. The temperature maps were acquired via the proton resonance frequency shift (PRFS) [44]. In particular, the calculation included masking of low SNR voxels, masking of the expected heating area and of areas exhibiting strong thermometry artifacts, baseline subtraction, and a second-order baseline drift correction. The MRI slices were positioned in four stacks, namely the focus, sagittal, near field, and far-field stacks. The focus stack contained three slices, all other stacks contained one slice. The focus stack was oriented perpendicular to the transducer axis, centered on the target point, and contained three slices to provide a complete view of the volumetric temperature distribution in the target ROI, with the center slice of the focus stack functioning as feedback for the MPC algorithm. The sagittal stack was oriented in sagittal orientation and was centered on the focus spot. It provides information on the vertical alignment of the focus. The near-field stack was positioned at the pig’s skin to monitor excess heat generation in the hydrogel spacer, indicative of excessive absorption. The far-field stack was positioned deep to the focal area.

Postinterventional imaging

After all hyperthermia treatments were complete, the perfusion of the pig’s hind legs was visualized. The purpose of this was to determine whether the hyperthermia treatments had unintentionally induced vascular shutdown and to localize the blood vessels in the vicinity of the treated target areas. This step was only performed after the hyperthermia treatments as the presence of contrast agents in the tissue of the animals would have interfered with MRI thermometry [68]. The employed technique was contrast-enhanced (CE) MRI, using the 3D T1-weighted FFE sequences described above for treatment planning. The images were acquired twice, both before and one minute after intravenous contrast agent injection (0.1 mmol/kg gadoteric acid, Dotagraf^®, Jenapharm GmbH & Co. KG, Jena, Germany). The parameters and field of view of the post-treatment scans were kept the same as the planning scans. The presence of non-perfused volumes (NPVs) was determined by subtracting pre- and post-injection CE MR images.

Controller architecture

The controller is designed to create a uniform temperature distribution in a user-defined target region of interest (target ROI) in the central slice of the MR-thermometry focus stack. Steering of the focus position is performed via electronic beam steering exclusively. This results in maximum diameter of 18 mm for the target area, which is used in this study. The controller follows the MPC paradigm and was described in detail in a previous article [66].

State-space model

The state-space model that will be used in the MPC scheme is based on the Pennes bioheat transfer equation (PBE) [69], incorporating heat loss via local perfusion, diffusion of heat and heating by an external source: (1) $$\begin{array}{c}{\rho }_t{C}_t\frac{\partial T}{\partial t}= \nabla ·\left(k\nabla T\right)-W_b{C}_b\left(T-T_a\right)+Q\end{array}$$(1) where ${\rho }_t$ is the density of the tissue [kg/m³], $k$ is the thermal conductivity [W/(°C·m)], $C_t$ and $C_b$ are the specific heat capacities of tissue and blood [J/(kg·°C)], $W_b$ is the perfusion rate [kg/(m³·s)], $T$ is the tissue’s temperature [°C], $T_a$ is the arterial temperature [°C], and $Q$ is the heating power density delivered to the tissue [W/m³]. It was assumed that ${\rho }_t,$ $C_t,$ and $k$ are homogeneous and constant.

The MPC regulates the temperature in a two-dimensional grid of voxels with grid spacing $h,$ that is, the in-plane voxel size, and the grid positions $(n,m).$ The diffusion term of the PBE was therefore discretized with the finite difference approximation of the Laplacian operator in two dimensions. (2) $$\begin{array}{c}\begin{array}{ccc}\nabla ·\left(k\nabla T\left(m,n\right)\right)=k\Delta T\left(m,n\right)& =& k\left(\frac{T\left(m+1,n\right)+T\left(m-1,n\right)+T\left(m,n+1\right)+T\left(m,n-1\right)}{h^2}\right.\\ & \\ \mathrm{ }& \mathrm{ }& -\frac{4·T\left(m,n\right)}{h^2})\end{array}\end{array}$$(2)

As the information on the tissue temperatures in adjacent parallel planes is time-shifted due to the sequential acquisition of temperature maps, out of plane diffusion is handled by introducing the apparent perfusion rate $W_a,$ which replaces $W_b.$

The sonication points that were made available to the controller were chosen to lie in the center of voxels inside the target ROI. These voxels constituted the target ROI’s edge and were arranged in a checker pattern everywhere else inside the target ROI. The heating power density $Q$ was discretized under the assumptions that energy dissipation is limited to the targeted voxel and that the absorption coefficient $\alpha$ [${\mathrm{m}}^3\mathrm{ }°\mathrm{C}\mathrm{ }{\mathrm{J}}^{-1}$] is homogeneous throughout the tissue. Attenuation of the ultrasound in the subcutaneous fat layer and in muscle was calculated with the attenuation coefficients ${\mu }_f$ and ${\mu }_m$ taken from [70] and the tissue layer thicknesses $d_f$ and $d_m$ measured along the transducer axis. With the sonication power allocated to grid position $(n,m),$ $p(n,m)$ [W], and the voxel volume $v$ [m³], $Q$ is described by (3) $$\begin{array}{c}Q\left(n,m\right)={\rho }_t·C_t·\exp \left(-\left({\mu }_f{d}_f+{\mu }_m{d}_m\right)\right)·\alpha ·\frac{p\left(n,m\right)}{v}\end{array}$$(3)

Finally, the temperature $T$ was replaced by the temperature elevation above $T_a,$ $T\mathrm{\prime }.$ The step size for the discretization in time is the acquisition time for a full set of thermometry slices, $\Delta t.$ At timestep $k+1,$ $T\mathrm{\prime }$ in the voxel at the grid position $(n,m)$ is therefore described by (4) $$\begin{array}{c}\begin{array}{ccc}{T}_{k+1}^{\mathrm{\prime }}\left(m,n\right)& =& A_0{T}_k^{\mathrm{\prime }}\left(m,n\right)+A_1\left(T_k^{\mathrm{\prime }}\left(m+1,n\right)\right.\\ & & +T_k^{\mathrm{\prime }}\left(m-1,n\right)+T_k^{\mathrm{\prime }}\left(m,n+1\right)\\ \mathrm{ }& \mathrm{ }& \left.+T_k^{\mathrm{\prime }}\left(m,n-1\right)\right)+q·p\left(n,m\right)\end{array}\end{array}$$(4) (5) $$\begin{array}{c}\begin{array}{ccc}{A}_0& =& \left(1-4·\frac{k}{{\rho }_t{C}_t{h}^2}\Delta t-\frac{W_a{C}_b}{{\rho }_t{C}_t}\Delta t\right)\\ A_1& =& \frac{k}{{\rho }_t{C}_t{h}^2}\Delta t\\ q& =& \frac{\exp \left(-\left({\mu }_f{d}_f+{\mu }_m{d}_m\right)\right)·\alpha ·\Delta t}{v}\end{array}\end{array}$$(5) where $A_0$ is the heat retention coefficient, which incorporates the loss due to perfusion and the diffusion out of the voxel, $A_1$ is the per-dynamic inter-voxel heat transfer coefficient, and $q$ is the target voxel heating coefficient [°C/W]. The sonication powers $p(n,m)$ are constrained to be greater than or equal to zero and to total no more than a specified maximum power $p_{\mathit{max}}$ at each timestep $k.$ The maximum power $p_{\mathit{max}}$ was set to 60 W in all experiments.

Cost function

The optimal distribution of heating power was determined via constrained optimization of the cost function. For this controller, it has been given the form (6) $$\begin{array}{c}{J}_t= \sum _{j=1}^H\sum _{\left(m,n\right)\in R}\frac{1}{{H·N}_R}{\left(\frac{T_{\mathit{obj}}^{\mathrm{\prime }}-T_{k+j}^{\mathrm{\prime }}\left(m,n\right)}{T_{\mathit{obj}}^{\mathrm{\prime }}}\right)}^2\end{array}$$(6) where $H$ is the control horizon, $\mathrm{R}$ is the set of grid positions belonging to the voxels inside the target ROI, $N_R$ is the number of voxels inside $R,$ and $T_{\mathit{obj}}^{\mathrm{\prime }}$ is the target temperature elevation. The control horizon determines how many future timesteps are considered in the optimization. For this study, $H$ was set to 5 to allow for heat exchange between the target area’s center and edge to be modeled while minimizing runtime.

Focal spot calibration and model identification

The transducer’s alignment with the MRI’s coordinate system was verified using multiple 16 s long test sonications at 60 W of acoustic output power before each experiment. During the test sonication, the HIFU focus was rapidly moved on a trajectory consisting of the transducer’s natural focus and eight points equally distributed on an 8 mm circle around the natural focus by electronic beam steering. Offsets from the intended centers of heating were corrected by adjusting the transducer’s position inside the HIFU table.

After alignment was verified, the MRI thermometry maps gathered during three additional test sonications were used to calibrate the state-space model’s parameters. The estimation of the $A_0$ and $A_1$ parameters are performed via linear regression, for which all voxels within a radius of 2 cm around the target area are included. The voxels’ temperatures at each timestep $k+1$ after transducer shutoff were used as the dependent variable. The respective voxel’s temperature at time $k$ (parameter $A_0$) as well as the sum of the nearest neighbors’ temperatures at time $k$ (parameter $A_1$) were used as the regressors. The model used, therefore, has the following form: (7) $$\begin{array}{c}{T}_{k+1}^{\mathrm{\prime }}\left(m,n\right)= A_0\left(T_k^{\mathrm{\prime }}\left(m,n\right)\right)\\ +A_1\left(T_k^{\mathrm{\prime }}\left(m+1,n\right)+T_k^{\mathrm{\prime }}\left(m-1,n\right)+T_k^{\mathrm{\prime }}\left(m,n+1\right)+T_k^{\mathrm{\prime }}\left(m,n-1\right)\right) \end{array}$$(7)

The linear regression was performed using the Python package scikit-learn, version 0.18.1 [71].

After the calculation of $A_0$ and $A_1$ is complete, the target voxel heating coefficient $q$ is calculated via linear regression. Here, the total temperature increase across all voxels in a radius of 20 mm around the natural focus at each timestep $k$ serves as the dependent variable and the retained fraction of the power added to the target slice is used as the regressor. This results in a linear model of the form (8) $$\begin{array}{c}\sum _{\left(m,n\right)\in M}{T}_k^{\mathrm{\prime }}(m,n)= q·P·\Delta t·\sum _{i=0}^k{(A_0+4A_1)}^k\end{array}$$(8) where $P$ is the acoustic output power selected for the sonication, $M$ is the set of grid positions of all voxels within a radius of 20 mm around the transducer’s natural focus and $(A_0+4A_1)$ is the fraction of power retained in the focus slice from one cycle to the next.

In the first two animals, $q$ was calculated by dividing the accumulated thermal load in the last temperature map of the test sonication by the total test sonication duration at the time of acquisition. However, the resulting absorption coefficients were notably higher than expected and exhibited large variability, prompting the development of the method described above.

Data analysis

Sonication metrics

The quality of the hyperthermic temperature control during the sonications was evaluated using parameters that were calculated from the temperature voxel values within the target ROI. The target ROI was placed in the center temperature mapping slice of the focus stack. The steady-state performance of the MPC was calculated from PRFS voxel temperature data acquired between 300 s and the end of the sonication. The used parameters were the mean temperature in the target ROI T_Avg the 10th and 90th percentile of the voxel temperature distribution $T_{90}$ and $T_{10},$ the standard deviation of the voxel temperatures in the target ROI ${\sigma }_{\hspace{0.17em}}$ and the cost finction J, $H$ of 1.

Spatial power allocation and equivalent voxel perfusion

The ability of the controller to counteract localized heat loss was investigated by superimposing the power applied to each voxel on the contrast-enhanced MRI scans that were acquired after the conclusion of the treatments. For this purpose, the power distribution of each treatment was averaged over the first minute after reaching the 300 s time point and the last minute before transducer shutoff. Additionally, the combined effect of heat diffusion and perfusion on the temperature of each voxel during cooldown was determined via linear regression. This equivalent voxel perfusion rate $W_e(n,m)$ was superimposed on the same MRI slices to verify the colocalization of local heat sinks and blood vessels. It was calculated by fitting the discretized PBE to the thermometry data acquired during the first minute after transducer shutoff while holding $k=0:$ (9) $$\begin{array}{c}{T}_{k+1}\left(n,m\right)=\left(1-\frac{W_e(n,m)C_b}{{\rho }_t{C}_t}\Delta t\right)T_k\left(n,m\right) \end{array}$$(9)

Analysis of the model error over time

The temperature distribution predicted by Equation (4)(4) $$\begin{array}{c}\begin{array}{ccc}{T}_{k+1}^{\mathrm{\prime }}\left(m,n\right)& =& A_0{T}_k^{\mathrm{\prime }}\left(m,n\right)+A_1\left(T_k^{\mathrm{\prime }}\left(m+1,n\right)\right.\\ & & +T_k^{\mathrm{\prime }}\left(m-1,n\right)+T_k^{\mathrm{\prime }}\left(m,n+1\right)\\ \mathrm{ }& \mathrm{ }& \left.+T_k^{\mathrm{\prime }}\left(m,n-1\right)\right)+q·p\left(n,m\right)\end{array}\end{array}$$(4) was calculated for each timestep from the PRFS temperature map and the applied heating pattern of the previous timestep. The resulting prediction was subtracted from the PRFS temperature map measured at the respective timestep to obtain the prediction error in each voxel at each timestep. The arithmetic mean of the prediction error in each voxel of the target ROI was calculated in each timestep to obtain the mean target ROI prediction error ${\mathrm{\Omega }}_{\mathit{ROI}}.$ Additionally, the prediction error in a sub-volume of each target ROI ${\mathrm{\Omega }}_{\mathit{sub}}$ was investigated. The sub-volumes consisted of the voxel in the target ROI to which most energy had been allocated throughout the sonication and those of its nearest neighbors which also contained an allowed sonication point. The behavior of ${\mathrm{\Omega }}_{\mathit{ROI}}$ and ${\mathrm{\Omega }}_{\mathit{sub}}$ was modeled as a linear combination of the prediction errors arising from a time-dependent perfusion parameter mismatch $\Delta w$ and an absorption parameter mismatch $\Delta q:$ (10) $$\begin{array}{c}\begin{array}{c}\mathrm{\Omega }\left(t,p\right)=\mathrm{ }\Delta w\left(t\right)+\Delta q·p_{\mathit{tot}}+{\mathrm{\Omega }}_0\\ \Delta w\left(t\right)=m_{w}t\end{array}\end{array}$$(10) where $p_{\mathit{tot}}$ is the total power applied to the respective area (either the entire ROI or the sub-ROI) and $m_w$ is the slope of the temporal change. The intercept ${\mathrm{\Omega }}_0$ incorporates the perfusion- and heat diffusion parameter errors arising from the a-priori model identification and the equilibration of diffusion losses. The diffusion losses were assumed to reach equilibrium by 300 s. Therefore, all the time points between 300 s and the end of the sonication were considered for this analysis. Experiments that were terminated prematurely and experiments that exhibited signs of phase drift were excluded.

Results

Model identification

The model parameters calculated from the PRFS temperature maps acquired during the test sonications and the derived thermal properties are listed in Table 1. The listed parameters were used as model parameters in the hyperthermia experiments following the respective model identification, see Section “Focal spot calibration and model identification”.

Table 1. Identified model parameters and corresponding physiological quantities.

Animal number/experiment number	$A_0$	$A_1$	$q$(°C/J)	$k$(W/m/°C)	$W_a$(kg/s/m³)	$\alpha$(mm³$·$°C/J)
1/1	0.520	0.120	[0.827]	0.39	0	[43.62]
2/1	0.500	0.110	[0.558]	0.35	14.51	[28.23]
2/2	0.510	0.110	[0.478]	0.35	12.09	[26.29]
3/1	0.513	0.117	0.397	0.38	4.37	21.27
3/2	0.500	0.117	0.330	0.38	7.64	17.43
4/1	0.540	0.110	0.403	0.35	4.84	20.46
4/2	0.570	0.100	0.356	0.32	7.25	17.76
4/3	0.503	0.117	0.374	0.38	7.01	19.53
5/1	0.507	0.107	0.295	0.35	15.70	16.34
5/2	0.533	0.110	0.335	0.35	6.53	16.87
Average (SD)	0.520 (0.021)	0.112 (0.006)	0.356 (0.036)	0.36 (0.02)	7.99 (4.58)	18.51 (1.76)

Note. Bracketed $q$ and $\alpha$ values were calculated using a previous algorithm utilizing only the final temperature map of the test sonication and are not included in the calculation of the respective averages. $A_0:$ heat retention parameter; $A_1:$ heat exchange parameter; $q:$ target voxel heating coefficient; $k:$ thermal conductivity; $W_a:$ apparent perfusion rate; $\alpha :$ absorption coefficient.

Hyperthermia experiments

The controller performance for all hyperthermia experiments is summarized in Table 2. Among $T_{\mathit{Avg}},$ $T_{10}$ and $T_{90},$ only $T_{90}$ showed a statistically significant difference between the value during the beginning of the steady-state phase ($t_{\mathit{Start}}$) and the value shortly before transducer shutdown ($t_{\mathit{Stop}}$) in a two-sided, paired t-test (p = .01). An example PRFS temperature curve and corresponding PRFS temperature maps for different phases of the hyperthermia sonication are shown in Figure 2. The controller first reached a mean temperature of 42 °C in the target area at 53.7 ± 9.0 s after the start of the sonication. The average absolute tracking error (${|T}_{\mathit{Avg}}-T_{\mathit{obj}}|$) during steady-state was 0.27 ± 0.15 °C. Treatment 4/3 was canceled prematurely due to time constraints in the animal protocol. NPVs were not observed in post-treatment CE MRI. Except for experiments 4/1 and 5/1, the average PRFS thermometry readings in the target ROI returned to baseline in all experiments, indicating the absence of excessive phase drift. The drift in experiments 4/1 and 5/1 was traced back to a problem in the phase unwrapping step.

Figure 2. Example of temperature evolution and temperature/power distribution during different control phases during a sonication, treatment 5/1. The black outline shows target ROI and the points inside the target ROI indicate power allocation. (A) Warmup phase, power is allocated primarily to the ROI’s center. (B) Beginning of steady-state phase; power is allocated primarily to the edges of the target ROI. (C) Late steady-state phase, immediately before the end of the sonication. Power is allocated to the edge and heterogeneously throughout the target ROI.

Table 2. Overview of sonication metrics.

Animal number/ experiment number	Animal weight [kg]	$p_{\mathit{out}}$(W)			$T_{\mathit{Avg}}$(°C)			$T_{90}$(°C)			$T_{10}$(°C)			$\sigma$(°C)
		Averaged timespan			Averaged timespan			Averaged timespan			Averaged timespan			Averaged timespan
		$t_{\mathit{tot}}$	$t_{\mathit{Start}}$	$t_{\mathit{Stop}}$	$t_{\mathit{tot}}$	$t_{\mathit{Start}}$	$t_{\mathit{Stop}}$	$t_{\mathit{tot}}$	$t_{\mathit{Start}}$	$t_{\mathit{Stop}}$	$t_{\mathit{tot}}$	$t_{\mathit{Start}}$	$t_{\mathit{Stop}}$	$t_{\mathit{tot}}$	$t_{\mathit{Start}}$	$t_{\mathit{Stop}}$	$J$$(\times 10^3$)
1/1	37	9.0	7.3	13.1	41.6	41.7	41.4	41.1	41.2	40.8	42.1	42.2	41.9	0.41	0.39	0.47	14
2/1	39	10.5	8.8	13.0	42.2	42.3	42.0	41.5	41.6	41.3	42.8	42.9	42.7	0.50	0.50	0.57	12
2/2		19.4	14.0	22.6	41.9	42.1	41.9	41.0	41.4	40.7	42.7	42.8	43.0	0.74	0.60	0.96	24
3/1	39	5.4	6.8	3.9	42.4	42.3	42.7	41.7	41.6	41.8	43.1	42.9	43.5	0.53	0.49	0.64	22
3/2		7.9	7.4	8.9	42.5	42.5	42.5	41.7	41.7	41.5	43.3	43.2	43.4	0.62	0.58	0.71	25
4/1	37	7.7	8.3	11.1	42.5	42.3	42.4	41.4	41.5	40.7	43.6	43.1	43.6	0.88	0.68	1.01	40
4/2		12.1	7.2	17.2	42.1	42.3	41.9	41.4	41.6	41.1	42.6	42.8	42.6	0.49	0.46	0.59	11
4/3		13.9	12.5	15.0	42.2	42.3	42.3	41.6	41.6	41.7	42.8	42.9	42.9	0.53	0.58	0.49	24
5/1	35	10.8	7.8	14.5	42.1	42.2	42.0	41.6	41.8	41.4	42.5	42.6	42.5	0.41	0.36	0.49	8
5/2		8.6	7.5	10.6	42.2	42.2	42.1	41.6	41.7	41.5	42.7	42.8	42.8	0.44	0.45	0.49	10
Average	37.4	10.5	8.8	13.0	42.16	42.19	41.97	41.46	41.58	41.26	42.82	42.80	42.87	0.557	0.517	0.650	19
SD	1.5	3.8	2.3	4.7	0.25	0.18	0.47	0.23	0.16	0.40	0.39	0.26	0.47	0.144	0.136	0.191	10

Note. The target temperature was always 42°C. p_out: HIFU output power; T_Avg: Mean temperature in target ROI; T₉₀, T₁₀: 90th and 10th percentile of voxel temperatures in the target ROI; σ: standard deviation of the voxel temperatures in the target ROI; t_tot: period from 300 s to transducer shutoff; t_Start: period from 300 s to 360 s; t_Stop: last minute before transducer shutoff; J: Cost function during t_tot.

Animals 2 and 3 exhibited localized erythema and swelling on the treated leg. The erythema was colocalized with air bubbles visible on MRI and resolved within 24 h in both animals. Limping or other behaviors indicating discomfort in the treated leg were not observed in any animal.

Spatial power allocation and equivalent perfusion

At the beginning of the sonications, heating power was allocated to the edges of the target ROI and was elevated in areas that showed elevated perfusion in the CE MRI scans (Figures 3 and 4, leftmost columns). By the end of the sonication, the latter effect was enhanced in almost all cases (Figures 3 and 4, middle columns). The only exception in this regard was treatment 4/3, which was canceled after only 10 min to stay below the maximum anesthesia duration defined in the animal protocol. The allocation of additional heating power at the end of the sonication was colocalized with areas of elevated $W_e(n.m),$ which were identified by fitting the thermometry data that was acquired after transducer shutoff (Figures 3 and 4, rightmost columns).

Figure 3. Overlay on contrast enhanced MRI images of the respective anatomical region, for animals 1, 2, and 3: Average power application pattern during the first minute after warmup (leftmost column), during the last minute before HIFU shutdown (middle column) and equivalent voxel perfusion during the first minute of cooldown (rightmost column).

Figure 4. Overlay on contrast enhanced MRI images of the respective anatomical region, for animals 4 and 5: Average power application pattern during the first minute after warmup (leftmost column), during the last minute before HIFU shutdown (middle column) and equivalent voxel perfusion during the first minute of cooldown (rightmost column).

Analysis of dynamic processes

The difference between time-dependent prediction errors $\Delta w$ of the sub-volumes and of the entire target ROIs between the equilibration at 300 s and the end of the sonication show a clear negative correlation with time, indicating a localized increase of heat loss over the period of the HIFU treatment (Figure 5(A)). Over the course of the 25 min between equilibration and the end of the sonication, five of seven power levels rose by a factor between 1.54 and 2.03. One power level rose by a factor of 2.81 (treatment 4/2) and one power level rose by a factor of 1.07 (treatment 3/2). On average, the power level required to maintain hyperthermia rose by a factor of 1.86 (Figure 5(B)). Changes in the mean absolute tracking error in the target ROI $\delta$ were small in five of the seven treatments with changes in the range of 0.1 K. Two treatments (1/1,2/2) showed an increase in $\delta$ of more than 0.2 K (Figure 5(C)). All sub-volumes except the sub-volume of treatment 4/2 contained a blood vessel with a diameter close to the in-plane voxel size of 1.53 mm. The closest vessel in treatment 4/2 was at a distance of 7.1 mm from the sub-volume and had a diameter of ∼2.6 mm, which is close to twice the in-plane voxel size (Figure 5(D–J)).

Figure 5. (A) The time-dependent model error in the sub-volumes around the voxel that had received the highest amount of heating energy minus the time-dependent model error in the entire ROI. (B) Total applied HIFU power normalized by the output power at 300 s. (C) Mean absolute tracking error in the target ROI. A one-minute moving average was applied to all data points in (A–C) to enhance clarity. (D–J) Contrast-enhanced MRI with superimposed outlines of the investigated sub-volumes in order from smallest to largest slope in (A). White arrowheads indicate blood vessels that adjoined or penetrated the target ROIs.

Discussion

We have investigated the performance of a model-predictive controller for MR-HIFU-mediated local hyperthermia at 42 °C for up to 30 min in a series of in vivo studies. The controller quickly stabilized the temperature in the target area with small steady-state tracking errors and narrow temperature distributions. While the heating power was uniformly distributed across the target area at the start of each experiment, the controller gradually allocated additional heating power to areas with local heat sinks, which were colocalized with blood vessels. Analysis of the difference between the modeling error at these localized heat sinks and the modeling error in the entire target ROI revealed a time-dependent component with a negative slope, which accumulated considerable magnitude in some cases.

In a previous article, we have shown that the MPC outperforms the state-of-the-art binary hyperthermia controller in a tissue-mimicking phantom both in terms of tracking error and uniformity of the temperature distribution [57,66]. A comparison of the in-vivo results presented here with the in-vivo results published for the binary controller that was acquired under similar conditions (target areas in the thigh muscle of pigs, 18 mm diameter, target temperature of 42 °C or 42.5 °C [57,72]) also indicates that the MPC generates a narrower temperature distribution in-vivo ($T_{10}-T_{90}=$ 1.25 ± 0.22 °C vs. 1.84 ± 0.23 °C for the MPC and the binary controller, respectively). In terms of target temperature tracking, the difference in the controllers’ performance is smaller than the measurement error (${|T}_{\mathit{Avg}}-T_{\mathit{obj}}|=\mathrm{ }$0.27 ± 0.14 °C vs. 0.19 ± 0.12 °C for the MPC and the binary controller, respectively) based on the data presented here and in the literature [57,72].

It is likely that the tracking errors of both controllers can be reduced by more closely adapting the applied heating power to the energy losses by diffusion and perfusion. In the MPC, this could be achieved by refining the fidelity of the model. The average thermal conductivity $k=0.36\pm 0.02\text{ W }{\mathrm{m}}^{-1}\mathrm{ }°{\mathrm{C}}^{-1}$ determined in the model calibration was lower than the value of $0.46\mathrm{W}\mathrm{ }{\mathrm{m}}^{-1}\mathrm{ }°{\mathrm{C}}^{-1}\mathrm{ }$ reported in the literature [73], likely resulting in an underestimation of in-plane diffusion. The average apparent perfusion parameter $W_a=7.99\pm 4.58\mathrm{kg}\mathrm{ }{\mathrm{s}}^{-1}{\mathrm{m}}^{-3}$ was much higher than the literature reference of $1.93\mathrm{kg}\mathrm{ }{\mathrm{s}}^{-1}{\mathrm{m}}^{-3}$ for regional blood flow [74]. While the latter was expected due to the incorporation of out-of-plane diffusion into $W_a,$ it may be possible to improve the target temperature tracking performance of the MPC slightly by directly modeling out-of-plane diffusion and implementing more refined parameter estimation methods such as the ones proposed by Dillon and Dragonu [75,76]. However, as hyperthermia is expected to stimulate perfusion, the identified perfusion parameter may become increasingly inaccurate [31,77]. The extension of the MPC to a fully offset-free MPC algorithm as recently developed by Deenen et al. [78] may therefore yield more robust results by constantly updating the model based on comparisons between the predicted and the measured temperature distributions of past feedback steps.

The MPC’s advantage over the binary algorithm in terms of uniformity may stem from the MPC’s flexible allocation of heating power, while the binary scheme only permits the application of one power level per feedback cycle throughout the target ROI. Consequently, the MPC algorithm can compensate for spatially heterogeneous absorption and heat sink effects more effectively and can react to the increase of perfusion that occurs during hyperthermia by adapting the used sonication pattern. Throughout this study, we have consistently observed an allocation of additional heating power to the locations of blood vessels, which became more pronounced over the course of the treatment and aligned with areas of rapid heat loss identified after transducer shutoff. This targeted compensation of local heat sinks should allow to adequately heat heterogeneously perfused tumors without causing vascular shutdown due to overheating, which is still a challenge today [19,79]. This is also made possible in part by the rapid feedback provided by PRFS thermometry, which allowed the controller to maintain highly uniform temperature distributions despite the pronounced local modeling errors that were observed in this study (Figure 5(A)).

Any results regarding perfusion obtained in this study must be considered qualitative due to the lack of a gold standard measurement like radioactive microspheres, and the mismatch of the calculated perfusion and diffusion parameters with values reported in the literature [73,74,80,81]. The CE MRI data in particular can only serve as an indicator for the presence of blood vessels as the vessel diameters are close to the voxel size of the images and the time periods between the hyperthermia treatments and the acquisition of the CE MR images were not standardized. Nonetheless, a thermoregulation response can clearly be seen via the localized change in the modeling error close to blood vessels relative to the rest of the ROI and the increase in the power required to maintain the target temperature. Working under the assumption that diffusion heat loss remains constant after 300 s, we can establish a lower bound for the thermoregulation response in the target area. In our study, we observed an increase in heating power by a factor of 1.86 between 300 s and transducer shutoff on average. Given that the power applied at 300 s was expended entirely to compensate for perfusion losses, this translates into an increase in perfusion by the same factor, which is close to the responses reported by Lokshina et al. for the leg muscles of male Sprague-Dawley rats [82]. As at least part of the power at 300 s must have been expended to compensate for diffusion losses, the true factor is likely higher than 1.86 in pigs. To our knowledge, comparable results have not yet been published for large animals or human subjects.

The visualization of spatially heterogenous heat sink effects and the estimation of perfusion changes as a response to hyperthermia presents an opportunity for future research. Most studies concerning the physiological response of tumor vasculature to elevated temperatures were conducted in the 1980s and only quantified averages across the whole tumor via radioactive microspheres [79]. Using a high-performance controller for MR-HIFU hyperthermia and post-hoc, high-fidelity thermal modeling on small length scales, it may be possible to gather new insight on the optimal conditions for the reperfusion and reoxygenation of tumor entities that are highly refractory to chemo- and radiotherapy. A suitable thermal model for the purpose of estimating the local ratio between perfusion and diffusion losses in heterogeneous tissues may be the discrete vasculature (DIVA) model originally developed by Kotte et al. [83,84]. In thermal simulations, the results obtained with models featuring discrete vasculature differ substantially from results obtained with models based on the Pennes bioheat equation [85] and the former show much better agreement with noninvasive thermometry data [86].

Future avenues for improvement include a more robust thermometry algorithm. This will become especially important for applications in abdominal organs, where both periodic and aperiodic motion will be encountered [87]. Possible strategies of implementing thermometry for MR-HIFU hyperthermia under motion have recently been demonstrated by Pichardo et al. in the form of a multi-baseline approach with the use of a pencil-beam navigator [88] and by Bing et al. in the form of thermometry acquisition during repeated breath holds [89]. For the heating of extended volumes under motion, the application of heat in a fractionated manner as recently demonstrated by Santos et al. appears to be the most feasible option in the short term [90] and appears to be synergetic with an MPC algorithm that incorporates sonication breaks in its state-space model. Furthermore, this control algorithm could be extended to treat larger target volumes by including temperature control in additional parallel monitoring slices and adding lateral transducer movement. With these additions and to safely treat highly perfused volumes, a strategy for near-field temperature management would be necessary as well [57]. Integration of these features will be indispensable to the clinical translation of the controller. Finally, we have observed considerable localized deviations from the expected thermal behavior of the heated region that is consistent with a thermoregulation response (Figure 5(A)). While this does not necessarily lead to a substantial decline in controller performance thanks to the rapid thermometry feedback available with MR-HIFU, it warrants the investigation of offset-free control strategies [78]. Besides improving controller performance both in terms of target temperature tracking and uniformity of the temperature distribution, these may be of special interest for drug delivery applications with temperature-sensitive liposomes, where the rate of release depends on the temperature of the blood [91].

The main limitation of this study is the lack of a tumor in the target area. The perfusion of tumors is heterogeneous [79], generally reacts less strongly to heating than muscle tissue [92,93] and considerable variability among preclinical tumor entities has been observed [94,95]. Tumors may also have acoustic characteristics distinct from healthy surrounding tissue [96]. The controller performance reported here may therefore not be representative of the performance that would be observed in the clinical case. Nonetheless, the developed MPC exhibits high performance in the presence of spatially heterogeneous and time-dependent perfusion. We, therefore, believe that the investigation of MPC algorithms for MR-HIFU hyperthermia in the clinical setting is warranted.

In conclusion, we have demonstrated the use of an MPC algorithm for MR-HIFU-mediated hyperthermia in a large animal model. The MPC algorithm provides uniform temperature distributions in the presence of heterogeneous perfusion by allocating additional heating power to localized heat sinks. MR-HIFU-mediated highly uniform hyperthermia holds the promise of detailed insight into the interplay between spatiotemporally heterogeneous perfusion and heating and the thermal enhancement effect. Clinically, it may become a tool for the treatment of deep-seated tumors that cannot be heated uniformly with other modalities.

Acknowledgments

The authors would like to thank Simon Hubertus (Profound Medical, Mississauga, Canada) for the acoustic field simulations of the test sonication protocol.

Disclosure statement

Edwin Heijman is employed by Philips. Sin Yuin Yeo is partly employed by Profound Medical. Juan Daniel Castillo Gómez is employed by Pentixapharm GmbH.

Additional information

Funding

This work was supported by the European Union FP7 Health program Health (“IPaCT”, grant agreement no. 603028) and the German Federal Ministry of Education and Research (“TSL-LIFU”, FKZ: 13XP5014C and “MR-HIFU-Pancreas”, FKZ: 13GW0364D).