Model predictive control for MR-HIFU-mediated, uniform hyperthermia

Abstract

Purpose: In local hyperthermia, precise temperature control throughout the entire target region is key for swift, safe, and effective treatment. In this article, we present a model predictive control (MPC) algorithm providing voxel-level temperature control in magnetic resonance-guided high intensity focused ultrasound (MR-HIFU) and assess the improvement in performance it provides over the current state of the art.

Materials and methods: The influence of model detail on the prediction quality and runtime of the controller is evaluated and a tissue mimicking phantom is characterized using the resulting model. Next, potential problems arising from modeling errors are evaluated in silico and in the characterized phantom. Finally, the controller’s performance is compared to the current state-of-the-art hyperthermia controller in side-by-side experiments.

Results: Modeling diffusion by heat exchange between four neighboring voxels achieves high predictive performance and results in runtimes suited for real-time control. Erroneous model parameters deteriorate the MPC’s performance. Using models derived from thermometry data acquired during low powered test sonications, however, high control performance is achieved. In a direct comparison with the state-of-the-art hyperthermia controller, the MPC produces smaller tracking errors and tighter temperature distributions, both in a homogeneous target and near a localized heat sink.

Conclusion: Using thermal models deduced from low-powered test sonications, the proposed MPC algorithm provides good performance in phantoms. In direct comparison to the current state-of-the-art hyperthermia controller, MPC performs better due to the more finely tuned heating patterns and therefore constitutes an important step toward stable, uniform hyperthermia.

Keywords:

• High intensity focused ultrasound
• noninvasive thermometry
• modeling
• control systems engineering/treatment optimization
• ;
• hyperthermia
• phantom
• tissue-mimick
• perfusion

Introduction

Mild local heating of tumor tissue to hyperthermic temperature in the range of 41–43 °C has been shown to have strong synergistic effects on radio- and chemotherapy, acting on cellular and tissue level [1–6]. Besides a plethora of preclinical evidence, several clinical trials demonstrated the synergistic effects of hyperthermia on radio-, chemo-, and chemoradiotherapy leading to improved outcome, such as local control, disease free-, and long-term survival [7–16]. As clinical studies also revealed that the outcome is sensitive to the applied thermal dose throughout the treatment, various measures have been introduced to describe and characterize the obtained temperatures over time and across the target tissue [17]. Induction of uniform hyperthermia across a treatment area is impeded by two major obstacles: first, the obtained temperatures depend on spatially and temporally varying tissue properties such as perfusion and heat diffusion. The former often leads to lower than intended heating close to blood vessels, while hyperthermia devices utilizing radiofrequency heating intrinsically lack the ability to compensate for local heat loss on the associated millimeter scale [18]. Second, temperature monitoring is currently performed using thermocouples inside or near the heated tissue [19], which only provides sparse spatial sampling of the temperature. Radiofrequency-based hyperthermia devices integrated into a magnetic resonance imaging (MRI) scanner could address some of above shortcomings and are currently evaluated for use in clinical settings [20,21].

Focused ultrasound was already explored in the 1930s and was employed early on for local noninvasive heating of deep-seated tissue to induce brain lesions [22,23]. Clinical application of high intensity focused ultrasound (HIFU), however, was only achieved after image guidance was added for spatial targeting and therapy control using diagnostic ultrasound or MRI [24–27]. The integration of an extracorporeal HIFU transducer into the patient bed of an MRI scanner (MR-HIFU) in particular enables both the in situ localization of the target tissue and the optimization of the transducer position before the treatment. Furthermore, the MRI provides near real-time temperature maps while the HIFU heating is ongoing. This is exploited for closed-loop feedback to the HIFU transducer, allowing to deliver a well-defined thermal dose and to maintain stable temperatures [28]. HIFU was recognized and tested as a hyperthermia device early on [24,29–31], but only MR-feedback enabled a stable hyperthermia over an extended period of time [32–34]. With a focus point of 2–3 mm in diameter and 5–7 mm in height, HIFU enables the local heating of tissue on the millimeter scale. For heating of larger volumes, electronic beam steering by coordinated ultrasound phase shifts was implemented to move the focus point over predefined circular perimeters [35]. In combination with MR thermometry and a simple binary feedback loop, electronic beam steering has been used to achieve MR-HIFU mediated, regional hyperthermia. The feasibility of this technique has been demonstrated in several preclinical studies [36–42] and recently in a first clinical trial using the MR-HIFU system Sonalleve^® (Profound Medical Inc., Toronto, Canada) [43]. By adding mechanical transducer movement facilitated by a robotic positioning system, a first large volume hyperthermia application based on MR-HIFU has been demonstrated in vivo [44] and applied clinically [45]. The feedback algorithm used in the referenced demonstration periodically repositions the transducer to heat different parts of the target region of interest (treatment cell) based on the lowest average temperature inside seven predefined sub-cells. Within the sub-cells, heating is regulated by adjusting the power delivered to the sub-cell as a whole and skipping points inside the sub-cell which are already above the target temperature. It is non-parametric in nature, i.e., it is designed to perform reliably without any prior knowledge of the particular tissue properties if they fall inside the anticipated range. Algorithms which do utilize tissue models for thermal therapies, however, can be more effective in achieving their respective control objective. Salomir et al. [32] as well as Quesson et al. [33], for instance, used physical models to calculate optimal HIFU intensities and thereby induced predefined temperature evolutions in the HIFU’s natural focus spot with high accuracy. Model predictive control (MPC) [46,47], a well-established control scheme using mathematical models of the to-be-controlled system in constrained optimization problems to calculate optimal control strategies, has been applied successfully as well: in ablation treatments, MPC has been shown to enable increased speed and patient safety in the delivery of a target thermal dose [48–50]. An MPC scheme for mild regional hyperthermia, which aims to directly control the temperature in a target ROI, has been introduced as well but was only demonstrated in silico thus far [51].

Here, we explore a novel control algorithm for MR-HIFU mediated mild hyperthermia applications, following the MPC paradigm. It aims to improve upon the existing binary control option by Tillander et al. [44], which we consider to be the state of the art, by optimizing the heat delivered to each individual voxel inside the target ROI based on a mathematical model of the tissue. As the appropriate model parameters for any individual case are generally unknown beforehand, the extent to which errors in the controller’s model can lead to problematic controller behavior was investigated in silico, followed by the quantification of the performance attainable with pretreatment model identification in phantom experiments. Since the controller was implemented on a clinical MR-HIFU system, a direct performance comparison between the MPC and the existing binary control algorithm of Tillander et al. [44] was also performed and is presented last.

Materials and methods

Model predictive control

MPC is a form of control that exploits a mathematical model of the system to be controlled (the state-space model) to predict its future state depending on the applied control actions. At each sampling instant, the state variables of the system, e.g., the temperatures in a set of voxels, are measured. This measurement is then used as a starting point to calculate the next control actions, e.g., a sequence of heating patterns, that minimize the objective function over a certain control horizon. The objective function’s minimum reflects the goals of the controller and depends on the predicted future state of the system. This calculation is usually done via constrained optimization, where the objective function is optimized by varying the control inputs while constraining the variables describing the system’s state (the state variables) to behave according to the model [46,47,52].

Controller architecture

State-space model

The state-space model that will be used in the MPC scheme is based on the Pennes bioheat transfer equation (PBE) [53] (1) $$\begin{array}{c}\rho C\frac{\partial T}{\partial t}= \nabla ·\left(K\nabla T\right)-W_b{C}_b\left(T-T_a\right)+\mathit{\text{Q}}\end{array}$$(1) where $\rho$ is the density of the tissue (kg/m³); $K$ is the diffusion coefficient (m²/s); $C$ and $C_b$ are the specific heat capacities of tissue and blood (J/(kg °C)), respectively; $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³).

Our implementation aims to control the temperature inside the target region using a single thermometry slice, which is acquired repeatedly alongside other slices that are required for monitoring purposes. A linear discrete-time state-space model of the temperature evolution in two spatial dimensions is therefore required. The discrete timestep size equals the acquisition time of one feedback cycle (i.e., the acquisition of a full set of thermometry slices) and the discrete space step size equals the voxel size in the acquired two-dimensional temperature maps. Summarizing diffusion and perfusion effects into dimensionless parameters, one obtains a simplified model for the predicted temperature elevation ${\widehat{T}}_{t+1}^{\prime }$ with (2) $$\begin{array}{c}{\widehat{T}}_{t+1}^{\mathrm{\prime }}\left(m,n\right)= A_0{T}_t^{\mathrm{\prime }}\left(m,n\right)\\ +A_1\left(T_t^{\mathrm{\prime }}\left(m+1,n\right)+T_t^{\mathrm{\prime }}\left(m-1,n\right)+T_t^{\mathrm{\prime }}\left(m,n+1\right)+T_t^{\mathrm{\prime }}\left(m,n-1\right)\right)\\ +\sum _{i=1}^{N_p}q·p_{t,i}·f_i(m,n)\end{array}$$(2) where $T_t^{\prime }(m,n)=T_t(m,n)-T_a$ is the temperature elevation above the arterial blood temperature $T_a$ at time $t\mathit{ }$ and discrete location $(m,n)$ (°C), $\mathit{ }{A}_0$ is the per-cycle heat loss coefficient, which incorporates the loss due to perfusion and the diffusion out of the voxel; and $A_1$ is the per-dynamic inter-voxel heat transfer coefficient. $q$ is the per-cycle, single-voxel temperature increase per Watt of acoustic power; $p_{t,i}$ is the acoustic power applied targeting the sonication point with index $i$ at timepoint $t;$ $f_i(m,n)$ is the distribution of heating power among voxel coordinates $\left(m,n\right)$ when targeting sonication point with index $i;$ and $N_p$ is the number of available sonication points. The discretization was performed under the assumptions that heat diffusion is homogenous in all observed voxels, the heat delivered in one feedback cycle is affected by diffusion and perfusion only from the next cycle on, and that heat is exchanged between nearest neighbors only.

The state variables are the temperature elevation $T_t^{\prime }(m,n)$ of the voxels at the locations $\left(m,n\right),(m,n)\in O$ at timepoint $t,$ where $O$ is the set of voxel coordinates inside the ROI which is passed to the controller (observed ROI). The input variables $p_{t,i},i=1,\dots ,N_p,$ are constrained to be greater than or equal to zero and to total no more than a certain maximum power $p_{\mathit{max}}$ at each timestep $t$ within the control horizon $H:$ (3) $$\begin{array}{c}0\le p_{t,i} \mathit{for} t=1,\dots ,H \mathit{and} i=1,\dots ,N_p\end{array}$$(3) (4) $$\begin{array}{c}\sum _{i=1}^{N_p}{p}_{t,i}\le p_{\mathit{max}}, \mathit{for} t=1,\dots ,H\end{array}$$(4)

Each allowed sonication point is chosen to coincide with the middle of an MRI voxel and is assumed to heat that particular voxel exclusively, i.e., $f_i(m,n)$ is assumed to be one at the targeted voxel’s coordinates and zero everywhere else.

Focal spot calibration and model identification

The transducer’s alignment with the MRI’s coordinate system was verified using multiple low-powered test sonications before each experiment. During the test sonication, the HIFU focus is moved by electronic beam steering on a trajectory consisting of the transducer’s natural focus point and eight points equally distributed on an 8 mm circle around the natural focus. Using these nine heated points as reference, errors in the transducer’s alignment are identified via spatial shifts from the expected heating pattern and are corrected by adjusting the transducer’s position inside the table accordingly.

The use of this multipoint test sonication pattern also creates a large number of warm voxels throughout the target area. The MRI thermometry maps gathered during the test sonications can therefore be used to calibrate the state-space model’s parameters to the target material. The estimation of the $A_0$ and $A_1$ parameters is 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: (5) $$\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}$$(5)

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

After the calculation of $A_0$ and $A_1$ is complete, the per-dynamic single-voxel heating rate $q$ is calculated by linear regression. Here, the total temperature increase across all voxels in a radius of 20 mm around the natural focus at each cycle $l$ 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: (6) $$\begin{array}{c}{\displaystyle \sum _{\left(m,n\right)\in M}}{T}_l^{\mathrm{\prime }}(m,n)= q·P·\sum _{i=0}^l{(A_0+4A_{\mathrm{i}})}^i\end{array}$$(6) where $P$ is the acoustic output power selected for the sonication, $M$ is the set of voxel coordinates 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.

Cost function

In MPC, the optimal control actions are determined via constrained optimization of the cost function, which is chosen to reflect the control objectives. For this controller, it has been given the form: (7) $$\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 }}-{\widehat{T}}_{t+j}^{\mathrm{\prime }}\left(m,n\right)}{5\mathbf{ }\mathrm{K}}\right)}^2\end{array}$$(7) where $H$ is the control horizon, $R$ is the set of coordinates belonging to the voxels inside the target ROI, $N_R$ is the number of voxels inside $R$, and $T_{\mathit{obj}}^{\prime }$ is the target temperature elevation.

Implementation

All phantom experiments were performed on a clinical MR-HIFU system (3T Achieva^®, Philips Healthcare, Best, The Netherlands, and Sonalleve^® V2 HIFU, Profound Medical, Toronto, Canada). We developed a custom MR-HIFU software providing treatment planning- and monitoring capabilities tailored to the MPC algorithm. The program was written in Python, version 2.7.13 [55]. PyQt 5.6.0 [56] was used to create the GUI and pyqtgraph 0.10.0 [57] provided efficient real-time visualization of the captured data. The communication with the MRI scanner and the HIFU system was achieved using the pyMRI and pyHIFU toolboxes [58]. The calculation of the optimal control actions was done using the python interface of Gurobi 7.5 [59]. The experiments employing the existing binary controller were performed using the Hyperthermia Release 3.5.955.1215 of the Sonalleve^® software (Profound Medical, Toronto, Canada).

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 = 296 mm, five slices, and voxel dimensions = 1.8 × 1.8 × 7 mm³, resulting in a dynamic scan time of 3.5 s (‘one cycle’) and a resolution of 160 × 160 voxels per slice. The temperature maps were acquired exploiting the proton resonance frequency shift (PRFS) [60]. The coils used for image acquisition were the HIFU table’s window coil and the HIFU Pelvis Coil (Model 905051-F, Philips Healthcare, Best, The Netherlands). In particular, the calculation included masking of low SNR voxels, masking of expected heating areas, 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 was positioned at the HIFU focus 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 controller. The sagittal slice was positioned to intersect the focus spot and provides information on the vertical alignment of the focus. The nearfield slice was positioned at the phantom’s acoustic window to monitor excess heat generation on the phantom’s surface or, in a hypothetical clinical setting, a patient’s skin. The far field slice was positioned near the ultrasound absorber at the far side of the phantom and is intended as a monitoring modality for organs and bone surfaces behind the focus.

Preparatory work

Phantom manufacturing process

All phantom experiments were performed with polyacrylamide hydrogels sealed in polypropylene containers. The containers can be opened on one side to provide an acoustic window and are lined with an ultrasound absorber on the opposite side to prevent reflections. The phantoms also contained a hollow channel of 5 mm diameter, which can be perfused with water to mimic a medium-sized blood vessel. The phantoms were manufactured in-house using a variation of the tissue-mimicking formulation by Negussie et al. [61]. The used chemicals were purchased from Sigma-Aldrich (Sigma-Aldrich Corporation, St. Louis, MO) and are detailed with their respective ratios in Table 1.

Table 1. Chemical composition of the phantoms used in this study.

Component	Product code	Fraction of volume
De-ionized water	–	82.24% (v/v)
Acrylamide/bis-acrylamide (ratio: 19:1), 40% (w/v) solution	A9926-5L	17.54% (v/v)
TEMED	T7024-100ML	0.21% (v/v)
APS	248614-500G	0.21% (w/v)
SiO2	S5631-1KG	2.00% (w/v)

TEMED: N,N,N0,N0-tetra-metyl-ethylenediamine; APS: ammonium peroxydisulfate; SiO₂: silicon dioxide particles, sizes 0.5–10 µm.

Model selection and identification

Ten test sonications with 60 W of acoustical power and 16 s duration were performed in the phantom to be used in the experiments, far away from the water channel. These test sonications employed the nine-point heating pattern described above. While the transducer was active and for 10 cycles after transducer shutoff, thermometry slices were acquired to provide training data for the thermal model.

As the detail of the inter-voxel heat exchange in the state-space model has a strong impact on the runtime of the controller, the chosen model must be simple enough to provide real-time control. The predictive quality of models with increasing numbers of nearest neighbors included in the prediction was therefore assessed and the peak runtimes required for solving the constrained optimization problem of the MPC at these different levels of detail was measured in simulations with 50 timesteps each. The simulations were run on the HIFU console (HP Z800 Workstation, Intel^® Xeon^® CPU X5650 @ 2.67 GHz) using the otherwise unchanged controller. The level of detail resulting in the shortest peak runtime while still providing markedly better predictive performance on the test set than all simpler models was selected for use in the remaining in silico and in vitro experiments. Half of all recorded voxel temperature measurements inside the observed 20 mm radius around the natural focus were assigned to the test set. Using the selected level of detail, the model parameters to be used in the following robustness experiments were determined.

Performance metrics

Controller performance was assessed in the steady state, which was established after 400 s at the latest in all experiments. The used metrics were the average steady-state performance $J,$ the average steady-state tracking error $T_{\mathit{Off}},$ and the average steady-state temperature distribution width $T_{\mathit{Range}}$ and are defined as (8) $$\begin{array}{c}J=\frac{1}{N_S·N_R}\sum _{t\in S}\sum _{\left(m,n\right)\in R}{\left(T_t\left(m,n\right)-T_{\mathit{obj}}\right)}^2\end{array}$$(8) (9) $$\begin{array}{c}{T}_{\mathit{Off}}=\frac{1}{N_S·N_R}\sum _{t\in S}\sum _{\left(m,n\right)\in R}{T}_t\left(m,n\right)-T_{\mathit{obj}}\end{array}$$(9) (10) $$\begin{array}{c}{T}_{\mathit{Range}}=\frac{1}{N_S}\sum _{t\in S}{T}_{90}\left(t\right)-T_{10}\left(t\right)\end{array}$$(10) where $S$ is the set of timepoints measured in the steady state ($400\mathit{ }\mathrm{s}<\mathit{ }t<600\mathit{ }\mathrm{s}$) and $N_S$ is the number of timepoints in $S.$ $T_{90}(t)$ and $T_{10}(t)$ denote the 90th and 10th percentile of voxel temperatures measured inside the target ROI at timepoint t.

Simulations

Influence of model errors in silico

In silico experiments allow the isolation of controller-specific errors from other sources of complications, such as alignment errors of HIFU- and MRI coordinate systems. The effect of plant-model mismatch, i.e., discrepancies between the model and the controlled system’s behavior, was therefore investigated in silico.

The controller’s model was initialized with a set of models containing varying $A_1$ and $q$ parameters while the properties of the simulated target tissue remained unchanged. The controller’s parameters were $q_{\mathit{MPC}}= (\raisebox{1ex}{$q_{\mathit{sim}}$}\!\left/ \!\raisebox{-1ex}{$3$}\right., q_{\mathit{sim}},3·q_{\mathit{sim}})$ and $A_{1,\mathit{MPC}}=(0.0,A_{1,\mathit{sim}},0.2)$ with respect to the simulated tissue’s parameters $q_{\mathit{sim}}$ and $A_{1,\mathit{sim}}.$ The simulations were performed by solving the finite difference Equation (2) on a two-dimensional grid of 1.8 mm spatial resolution and a 3.5 s temporal resolution. The used model was selected beforehand using thermometry data acquired during test-sonications on the polyacrylamide phantom described above. For these simulations, the controller was operated on a circular target ROI with a diameter of 18 mm. $T_{\mathit{obj}}$ was set to 5 °C above the baseline. To simulate noise, normally distributed random numbers with a standard deviation of 0.4 K were added to the temperature maps passed to the controller as feedback.

Phantom experiments

Influence of model errors in vitro

Nine hyperthermia experiments were performed on the previously characterized polyacrylamide phantom using MPCs initialized with all possible combinations of the average value found for $A_1,$ the case of negligible perfusion ($A_1=0.0)$, and immediate leveling of temperature gradients between nearest neighbors ($A_1=0.2)$ with the average value found for $q$ and $q$-mismatches by the factors $\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$ and 3. Heat loss other than in-plane diffusion was excluded from the starting model by setting $A_0=1 - 4A_1.$ Each hyperthermia experiment was run for 10 min with $T_{\mathit{obj}}=5 \mathrm{K},$ $P_{\mathit{max}}=60 \mathrm{W}$ and $H=5.$ For these experiments, the target area was placed far from the water channel and the channel’s perfusion was disabled.

Performance comparison: MPC and binary controller

The MPC’s performance was assessed by comparison with the current state-of-the-art controller introduced by Tillander et al. [44], which will be referred to as the binary controller. For the comparison in a homogeneous target area, the water supply of the phantom was disabled and a round target ROI with 18 mm diameter was placed several centimeters away from the water channel. For the comparison in an inhomogeneous target area, the phantom was perfused with de-ionized water at 2.4 ml/s and the target ROI was positioned 1 mm away from the water channel. Both target ROIs were sonicated with the MPC as well as the binary control algorithm. For both experiments, the MPC was initialized with parameters identified from test sonication data acquired before the experiment in the respective configuration. Each hyperthermia experiment was run for 10 min with $T_{\mathit{obj}}=5 \mathrm{K},$ $P_{\mathit{max}}=60 \mathrm{W}$ and $H=5.$

Results

Preparatory work

Model selection and identification

First, the coefficient of determination as well as the controller runtimes was measured using increasing numbers of nearest neighbors in the prediction of voxel temperatures (Figure 1). The model incorporating four nearest neighbors (N = 4) shows a significantly higher coefficient of determination than the model neglecting diffusion (N = 0) and led to a peak runtime of 390 ms. The model with eight nearest neighbors required a 55.6% longer peak runtime of 608 ms and improved the predictive quality only slightly ($R^2$=0.946 vs. 0.947 on average). Models with more nearest neighbors did not further improve the coefficient of determination significantly but only increased peak run time. The chosen model therefore incorporated four nearest neighbors only. The $A_1$ and $q$ parameters derived from the test sonications using this model were 0.126 ± 0.008 and 0.882 ± 0.094, respectively.

Figure 1. Coefficient of determination on the test set vs. the number of nearest neighbors used in the model and corresponding peak runtimes required for optimization.

Simulations

Influence of model errors in silico

The simulated temperature evolutions in the target ROI for plant-model mismatches in $A_1$ and $q$ are displayed in Figure 2. Regardless of the used $A_{1,\mathit{MPC}},$ the controller performs best in the cases with $q_{\mathit{MPC}}=q_{\mathit{sim}}.$ In the case of severely underestimated heating ($q_{\mathit{MPC}}=q_{\mathit{sim}}/3$), the temperature distributions are much broader than in the matched heating case for all $A_1$ parameters and show heat spikes throughout the duration of the sonication. Additionally, positive tracking errors between 304 and 542 mK are observed. In the overestimated heating case ($q_{\mathit{MPC}}=3·q_{\mathit{sim}}$), the widths of the temperature distributions are larger than for the matched heating case as well for all $A_1$ parameters and negative tracking errors between 285 and 353 mK are observed. The best performance $J$ was achieved using the correct model parameters.

Figure 2. Simulated temperature curves in the target ROI vs. time, not including the noise added to the feedback signal. In each simulation, the controller was given a different model while the simulated tissue remained unchanged. The eight panels around the center represent the mismatched cases with the used parameter displayed at the respective row and column. T_Range: average steady-state temperature distribution width; T_Off: average steady-state tracking error; J: average steady-state performance.

The controller’s responsiveness to measured temperature deficits is shown in Figure 3 for varying model parameters. Over- and underestimation of the $q$ parameter leads to under- and overcompensation of measured temperature deficits, respectively. For the case of neglected diffusion ($A_{1,\mathit{MPC}}=0.$0), the controller aggressively counteracts measured temperature deficits while in the case of immediate leveling of temperature gradients ($A_{1,\mathit{MPC}}=0.2$), the controller is much less responsive.

Figure 3. Temperature change of voxels inside the target ROI vs. measured temperature deficit at the previous timestep in silico. The eight panels around the center represent the mismatched cases with the used parameter displayed at the respective row and column. The used q parameter determines the magnitude of the response while the used A1 parameter influences the stringency with which measured temperature deficits are compensated.

Phantom experiments

Influence of model errors in vitro

The evolution of the temperature distributions during nine in vitro experiments with varying model parameters is shown in Figure 4. Across all experiments, the average absolute steady-state tracking error $T_{\mathit{Off}}$ was 138 mK and the average steady-state width of the temperature distribution $T_{\mathit{Range}}$ was 740 mK. In all tested cases, an underestimation of the heating power led to $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ values above 144 mK and 754 mK while an overestimation led to $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ values below –185 mK and above 785 mK. Experiments using the experimentally identified $q$ parameter performed better in both aspects with $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ values smaller than 62 mK and 676 mK. The use of the experimentally identified $A_1$ parameter led to the best performance $J$ for two out of three different $q$ parameters and the second-best result in one case. The best performance $J$ was achieved using the experimentally identified parameters.

Figure 4. Temperature curves in the target ROI vs. time for varying controller models in nine otherwise identical phantom experiments. The used values for the A₁ and Q parameters are given for each row and column. T_Range: average steady-state temperature distribution width; T_Off: average steady-state tracking error; J: average steady-state performance.

Performance comparison: MPC and binary controller

The MPC’s performance was rated in a direct comparison to the current state-of-the-art controller for MR-HIFU mediated mild hyperthermia. The comparison was made both for a homogenous target ROI (Figure 5) and a target ROI adjacent to the water channel inside the phantom material (Figure 6). The steady-state (400 s<t < 600 s) performance metrics $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ for the homogeneous case were 67 mK and 698 mK for the MPC and 605 mK and 1734 mK for the binary controller, respectively. The steady-state performance metrics $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ for the inhomogeneous case were 49 mK and 862 mK for the MPC and 646 mK and 1629 mK for the binary controller, respectively. The average spatial temperature profiles through the centers of the target ROI in steady state (Figures 5(B) and 6(B)) showed bell shapes for the binary controller and plateaus for the MPC controller. The MPC allocated heating power mainly to the edges of the target ROI and increased heating at the edge of the ROI bordering the water channel in the inhomogeneous case. Information on the power distribution used by the binary controller could not be obtained from the clinical software and is therefore not shown here.

Figure 5. (A) Temperature distribution vs. time achieved with the binary controller and the MPC algorithm in a homogeneous region of the phantom. (B) Spatial temperature profiles and the MPC’s power distribution profile through the center of the target region, average for 400 s<t < 600 s. Target area is shaded.

Figure 6. (A) Temperature distribution vs. time achieved with the binary controller and the MPC algorithm near the water channel of the phantom. (B) Spatial temperature profiles and the MPC’s power distribution profile through the center of the target region perpendicular to the water channel, average for 400 s<t < 600 s. Target area (grey) and water channel (purple) are shaded.

Discussion

The impact of errors in the estimation of tissue properties was explored in silico and in vitro. Both investigations showed that mismatches of the $A_1$ parameter have a weaker influence on steady-state controller performance compared to the influence of the heating parameter $q.$ This observation can be explained by the diminished influence of diffusion inside the target ROI once a highly homogeneous temperature distribution has been reached. Neglect of diffusion in the model ($A_1=0)$ together with a severe under- or overestimation of $q,$ however, led to increased tracking errors of up to 542 mK in silico and 230 mK in vitro, caused by an aggressive overcompensation of measured temperature deficits. While these errors are smaller than the tracking errors that were encountered using the binary controller, this still shows the necessity of modeling diffusion for precise temperature control. At the same time, modeling diffusion between a large number of voxels offers only small benefits beyond the already very good predictive performance achievable with four nearest neighbors.

Besides the expected tracking errors, the underestimation of $q$ caused heat spikes which were observed prominently in silico, where the feedback noise could be masked. This influence of the mismatch in $q$ arises from the resulting controllers’ overcompensation of measured temperature deficits induced by noise. Higher heating estimates, on the other hand, will lead to insufficient heating in response to any measured temperature deficit. This led to a smoother temperature evolution in the corresponding simulations, but not to an improvement in control performance over the matched case. One avenue that could be pursued in the future to reduce tracking errors and the controller’s responsiveness to noise in a controlled manner is the extension of the MPC algorithm to a full offset-free implementation, which is used in other control problems where the system’s state cannot be measured directly or the controller input is noisy [52,62].

In the direct comparison between the MPC and the binary control algorithm, the MPC implementation outperformed the binary controller both in terms of the tracking error and width of the temperature distribution after establishment of the steady state. In the homogeneous case, the MPC’s $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ were only 11.1% and 40.3% of the binary controller’s $T_{\mathit{Off}}$ and $T_{\mathit{Range}}.$ The reason for this becomes evident in the spatial temperature profiles produced by the two controllers (Figure 5(B)). The MPC’s profile shows a flat plateau inside the target area, while the binary controller gives rise to a bell shape, indicative of excess heating in the center of the target ROI. The MPC avoids this problem by optimizing the applied heating pattern at each time step, thereby shifting all applied power to the edge of the controlled area automatically after the ROI center is sufficiently heated. In the inhomogeneous case, where the heated region is placed adjacent to a perfused water channel, both controllers produce a skewed distribution (Figure 6(B)). The MPC’s advantage is considerable in this case as well with $T_{\mathit{Off}}$ and $T_{\mathit{Range}}$ at 7.6% and 52.9% of the binary controller’s $T_{\mathit{Off}}$ and $T_{\mathit{Range}},$ respectively. Analogous to the homogeneous case, this is a consequence of the flexible heating patterns used by the MPC, allowing to allocate additional heating power to the edge of the target ROI bordering the water channel. This capability might prove particularly useful in a clinical setting where regions of enhanced perfusion or variations in ultrasound absorption might not be detectable in advance or might arise during the treatment, e.g., by the heat-induced stimulation of blood flow.

The controlled environment that was used to show the MPC’s advantages over the current state-of-the-art in the absence of such unpredictable effects and to assess the problems particular to this new MPC algorithm is at the same time a limitation of this study. It cannot demonstrate adequately how the MPC would perform when encountering problems that are particular to the in vivo setting. Most of these problems, first and foremost those related to patient safety, e.g., protection of the skin, the detection of thermometry artifacts and detection of patient motion, must be addressed by appropriate patient preparation, patient monitoring and dedicated safety software modules separate from the controller, however. Nonetheless, validation of the controller’s performance in an in vivo setting, where perfusion might change over time or tissue properties may vary spatially in unpredictable ways, is required and will be performed in a dedicated study.

In conclusion, we presented the structure, implementation and in vitro evaluation of an MPC algorithm for improved uniform hyperthermia using MR-HIFU, which employs personalized thermal models derived from low-powered test sonications for each treatment. The MPC outperforms the current state-of-the-art controller in vitro both in terms of target tracking performance as well as spatial homogeneity of the resulting temperature distribution. This is made possible via the automatic and flexible redistribution of heating power by repeated model-based optimization. This also holds in the presence of a localized heat sink, which was counteracted by the MPC via automatic allocation of additional heating power. We therefore believe that the proposed MPC solution will further advance MR-HIFU based hyperthermia and may be key for applications where tight control of the temperature is essential, such as device-targeted drug delivery.

Disclosure statement

Edwin Heijman is employed by Philips.

Additional information

Funding

This study is supported by the European Union FP7 Health Program Health [“IPaCT”, Grant Agreement No. 603028], the German Federal Ministry of Education and Research [“TSL-LIFU”, FKZ: 13XP5014C], the Dutch Cancer Society and the Netherlands Organisation for Scientific Research (NWO) as part of their Joint Partnership Programme: “Technology for Oncology”. This project is partially financed by the PPP Allowance made available by Top Sector Life Sciences & Health. We thank Samuel Pichardo for granting access to the pyMRI and pyHIFU toolboxes.