- Email: [email protected]
Offset-free model predictive control for enhancing MR-HIFU hyperthermia in cancer treatment⁎
6th IFAC Conference on Nonlinear Model Predictive Control 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, USA, August 19-22, 2018 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, USA, August 19-22, 2018 Available online at www.sciencedirect.com Madison, WI, USA, August 19-22, 2018 6th IFAC Conference on Nonlinear Model Predictive Control Madison, WI, USA, August 19-22, 2018
ScienceDirect
IFAC PapersOnLine 51-20 (2018) 191–196
Offset-free model predictive control for Offset-free model predictive control for Offset-free model predictive control for enhancing MR-HIFU hyperthermia in Offset-free model predictive control for enhancing MR-HIFU hyperthermia in enhancingcancer MR-HIFU hyperthermia in ⋆ treatment ⋆⋆ enhancingcancer MR-HIFU hyperthermia in treatment cancer treatment ⋆ cancer ∗ ∗∗ ∗ D.A. Deenen ∗∗ E. Maljaarstreatment ∗ L. Sebeke ∗∗ Bram de Jager ∗
D.A. Deenen ∗ E.∗∗∗ Maljaars ∗ L. ∗∗ Sebeke ∗∗ Bram de Jager ∗ ∗ Heijman H. Gr¨ ullL. W.P.M.H.Bram Heemels D.A. E. Deenen E.∗∗∗ Maljaars de Jager ∗∗ Sebeke ∗ E. Heijman H. Gr¨ u Heemels ∗ ∗ll ∗∗ W.P.M.H. ∗∗ ∗ ∗∗∗ ∗ D.A. E. Deenen E. Maljaars Bram de Jager Heijman H. Gr¨ ullL. Sebeke W.P.M.H. Heemels ∗∗∗ ∗∗ ∗ ∗ E. Heijman H. Gr¨ ull W.P.M.H. Heemels ∗ Eindhoven University of Technology, Mechanical Engineering, University of Technology, Mechanical Engineering, ∗ Eindhoven Control Systems Technology, Eindhoven, The Netherlands Eindhoven University of Technology, Mechanical Engineering, ∗ Control Systems Technology, Eindhoven, The Netherlands Eindhoven University of [email protected]) Technology, Mechanical Engineering, (e-mail: Control Systems Technology, Eindhoven, The Netherlands (e-mail: [email protected]) ∗∗ Control Systems Technology, Eindhoven, TheGermany Netherlands (e-mail: [email protected]) Hospital Cologne, Cologne, ∗∗ University University Hospital Cologne, Cologne, Germany ∗∗ ∗∗∗ [email protected]) Philips (e-mail: Research Germany, University Hospital Cologne, Aachen, Cologne, Germany Germany ∗∗∗ Philips Research Germany, ∗∗ ∗∗∗ University HospitalGermany, Cologne, Aachen, Cologne, Germany Germany Philips Research Aachen, Germany ∗∗∗ Philips Research Germany, Aachen, Germany Abstract: Abstract: In In this this paper, paper, an an offset-free offset-free MPC MPC scheme scheme is is developed developed to to improve improve performance performance and robustness of temperature control inMPC magnetic focused Abstract: In this paper, an offset-free schemeresonance-guided is developed to high-intensity improve performance and robustness of temperature control in magnetic resonance-guided high-intensity focused Abstract: In this paper,cancer an offset-free is design developed to high-intensity improve and robustness of temperature control inMPC magnetic resonance-guided focused ultrasound hyperthermia treatments. The scheme proposed effectively handlesperformance constraints ultrasound hyperthermia cancer treatments. The proposed design effectively handles constraints robustness of temperature control to in achieve magnetic resonance-guided focused and exploits the body’s thermal faster and moreeffectively evenlyhigh-intensity distributed heating of ultrasound hyperthermia cancerbehavior treatments. The proposed design handles constraints heating of and exploitshyperthermia the body’s thermal behavior to achieve faster and moreeffectively evenly distributed ultrasound cancertreatment treatments. The proposed design handles constraints the tumor leading to improved with respect to the methods currently applied and exploits the body’s thermal behavior toquality achieve faster and more evenly distributed heating of the tumor leading to improved treatment with respect to the methods currently applied and exploits the proposed body’s thermal behavior toquality achieve faster andestimation more evenly heating of in The scheme incorporates disturbance todistributed remove the steadytheclinics. tumor leading to improved treatment quality with respect to the methods currently applied in clinics. The proposed scheme incorporates disturbance estimation to remove the steadythe tumor to from improved treatment qualitydisturbance with respect to given the methods currently in clinics. The proposed scheme incorporates estimation to remove the applied steadystate offsetleading resulting plant-model mismatch, which for the healthcare application is state offsetThe resulting from scheme plant-model mismatch, which for the given healthcare is in clinics. proposed incorporates disturbance estimation to removeapplication the steadyinevitable practice due the diversity and Simulation state offsetin mismatch, which for in thepatients given healthcare application is inevitable inresulting practice from due to toplant-model the large large physiological physiological diversity in patients and tumors. tumors. Simulation state offset resulting from mismatch, which for in thepatients healthcare application is results are in presented demonstrate algorithm’s capability togiven improve performance. inevitable practice to due toplant-model the largethis physiological diversity androbust tumors. Simulation results are presented demonstrate algorithm’s capability to improve performance. inevitable practice to due to the largethis physiological diversity in patients androbust tumors. Simulation results are in presented to demonstrate this algorithm’s capability to improve robust performance. © 2018,are IFAC (International Federation ofthis Automatic Control) Hosting to by improve Elsevier Ltd. All rights reserved. results presented to demonstrate algorithm’s capability robust performance. Keywords: Keywords: Healthcare, Healthcare, hyperthermia, hyperthermia, offset-free offset-free MPC, MPC, cancer cancer treatment, treatment, high-intensity high-intensity Keywords: Healthcare, hyperthermia, offset-free MPC, cancer treatment, high-intensity focused ultrasound. focused ultrasound. Keywords: Healthcare, hyperthermia, offset-free MPC, cancer treatment, high-intensity focused ultrasound. focused 1. ultrasound. INTRODUCTION via volumetric temperature maps obtained using an MRI 1. via volumetric temperature maps using an 1. INTRODUCTION INTRODUCTION via volumetric temperature maps obtained obtained usingallows an MRI MRI scanner, see Maloney and Hwang (2015). This for scanner, see Maloney and Hwang (2015). This allows for 1. INTRODUCTION via volumetric temperature maps obtained using an MRI a completely noninvasive procedure, thereby significantly scanner, see Maloney and Hwang (2015). This allows for Mild local hyperthermia involves the heating of a specific a completely noninvasive procedure, thereby significantly Mild local hyperthermia involves the heating of a specific see Maloney Hwang (2015). This allows the for to patientand well-being. However, realizing acontributing completely noninvasive procedure, thereby significantly Mild local hyperthermia the heating of a specific target volume inside the involves human body to temperatures of scanner, contributing to patient well-being. However, the target inside human body to temperatures of ◦volume completely noninvasive procedure, thereby realizing significantly desired temperature distribution and accurately maintainto patient well-being. However, realizing the Mild local hyperthermia the heating of a specific target inside the the involves human body toNumerous temperatures of acontributing C for approximately 90 minutes. clinical 39-45 ◦volume temperature distribution and accurately maintain39-45 ◦ C for approximately 90 minutes. Numerous clinical desired to course patient However, realizing the desired temperature distribution and accurately ing it over the ofwell-being. an entire treatment in maintaina clinical targethave volume inside human bodyisto temperatures of contributing 39-45 trials thatthe hyperthermia aNumerous potent adjuvant C forshown approximately 90 minutes. clinical ing it over the course of an entire treatment in a clinical trials have shown that hyperthermia is a potent adjuvant ◦ desired temperature distribution and accurately maintainsetting is no simple task. Combined with the fact that ing it over the course of an entire treatment in a clinical 39-45 C for approximately 90 minutes. Numerous clinical therapy for shown cancerthat treatment, see for example Datta et al. setting is no simple task. Combined with the fact that trials have hyperthermia is a potent adjuvant therapy for shown cancerthat treatment, see for example Datta et al. ing over thesimple coursetask. entire therapies, treatment in fact a of clinical the itbeneficial effects of an thermal and mild is no Combined with the that trials hyperthermia isal.a (2018). potent adjuvant (2015);have Mallory et al. (2016); Issels It signiftherapy for cancer treatment, see foretexample Datta et al. setting the beneficial effects task. of thermal therapies, andfact of mild (2015); Mallory et al. (2016); Issels et al. (2018). It signifsetting is no in simple Combined with the that hyperthermia particular, are closely related to the the beneficial effects of thermal therapies, and oftissue mild therapy for cancer treatment, seeprimary foretexample Datta et al. icantly increases the efficacy of methods such as (2015); Mallory et al. (2016); Issels al. (2018). It signifin particular, are closely related to the icantly increases the efficacy of primary methods such as hyperthermia beneficialtruly of thermal therapies, oftissue mild hyperthermia ineffects particular, are closely related and tosee the tissue temperatures achieved during treatment, Franck(2015); Mallory etthe al.efficacy (2016); Issels et al. (2018). signifchemo- increases and radiotherapy by sensitizing the treatedItsuch tumoricantly of primary methods as the temperatures truly achievedare during treatment, Franckchemoand radiotherapy by sensitizing the treated tumorin particular, closely related tosee the tissue ena this strongly motivates the temperatures truly achieved during treatment, see Franckicantly increases efficacy of methods as hyperthermia chemoandAlso, radiotherapy by the treated such tumorous tissue. itthe enables thesensitizing useprimary of temperature-sensitive ena et et al. al. (2009), (2009), this strongly motivates the development development ous tissue. Also, it enables the use of temperature-sensitive trulythis achieved during treatment, see Franckof feedback controllers for MR-HIFU. ena et al. (2009), strongly motivates the development chemoandAlso, radiotherapy the treated tumorliposomes for heat-mediated drug release, allowing for temperatures ous tissue. it enablesby thesensitizing use of temperature-sensitive feedback controllers for MR-HIFU. liposomes for heat-mediated drug release, allowing for of et al. (2009), this strongly motivates the development of feedback controllers for MR-HIFU. ous Also, it enables the anticancer use of temperature-sensitive moretissue. efficient of the drugsallowing specifically liposomes for delivery heat-mediated drug release, for ena Currently, (experimental) clinical more efficient delivery of the anticancer drugs specifically feedback (experimental) controllers for MR-HIFU. Currently, clinical applications applications of of local local hyhyliposomes forregion heat-mediated drug release, allowing for of to theefficient tumor at lower drug concentrations, more delivery of thesystemic anticancer drugs specifically perthermia(experimental) by MR-HIFU typically involve a (mostly) Currently, clinical applications of localprehyto the tumor region at lower systemic drug concentrations, perthermia by MR-HIFU typically involve a (mostly) premore efficient delivery of the anticancer drugs specifically to the tumor region at lower systemic drug concentrations, as in Hijnen et al. (2014); Staruch et al. (2015). Most Currently, clinical applications of local hydetermined(experimental) sonication plan, possibly extended with simple perthermia by MR-HIFU typically involve a (mostly) preas in Hijnen et al. (2014); Staruch et al. (2015). Most sonication plan, possiblyinvolve extended with simple to the tumor the region lower drug importantly, showsystemic that local does determined as in Hijnen et results al. at(2014); Staruch ethyperthermia al.concentrations, (2015). Most perthermia by MR-HIFU typically a (mostly) preonline MR-driven feedback algorithms. Examples include determined sonication plan, possibly extended with simple importantly, show Staruch that localethyperthermia does online MR-driven feedback algorithms. Examples include as Hijnenthe et results al. (2014); al. (2015). Most notinincrease toxicity inshow healthy tissuehyperthermia and the thereby importantly, the results that local does determined sonication plan, algorithms. possibly extended with and/or simple binary MR-driven strategies scale power online feedback Examples not increase toxicity in healthy tissue and the thereby strategies to to scale the the sonication sonication power include and/or importantly, results that local hyperthermia does induced side the effects, making it highly clin- binary not increase toxicity inshow healthy tissueappealing and the for thereby online feedback algorithms. Examples include length MR-driven of the heating interval (Enholm et al. (2010); Parbinary strategies to scale the sonication power and/or induced side effects, making it highly appealing for clinof the heating interval et (2010); Parnot toxicity in healthy tissueappealing and therates thereby induced side effects, making it highly for clinical increase application to increase treatment success and length binary strategies to basic scale the(Enholm sonication power and/or length of the(2012)), heating interval (Enholm et al. al. (2010); Partanen et al. proportional-integral-derivative ical application to increase treatment success rates and tanen et al. (2012)), basic proportional-integral-derivative induced side effects, making itpatients highly in appealing for clinimprove quality of life for the therapy. ical application to increase treatment success rates and length ofal. the(2012)), heatingbasic interval (Enholm al. (2010); Par(PID)-based methods (Mougenot et al.et(2009)) or some tanen et proportional-integral-derivative improve quality of life for the patients in therapy. methodsbasic (Mougenot et al. (2009)) or some ical application to life increase rates and (PID)-based improve quality of for thetreatment patients insuccess therapy. tanen al. (2012)), proportional-integral-derivative hybridetform ofmethods PID and(Mougenot bang-bang (Staruch al. (PID)-based etcontrol al. (2009)) or et some A particularly interesting technology for hyperthermia is form ofmethods PID and(Mougenot bang-bangetcontrol (Staruch et al. improve qualityinteresting of life for the patientsfor in therapy. A particularly technology hyperthermia is hybrid (PID)-based (2009)) some (2011); form Lin of et PID al. (1990); Bing et control al.al.(2015)). Aor major hybrid and bang-bang (Staruch et al. magnetic-resonance-guided high-intensity focused ultraA particularly interesting technology for hyperthermia is Lin et al. (1990); Bing et al. (2015)). A major magnetic-resonance-guided high-intensity focused ultra- (2011); form of and bang-bang (Staruch ettake al. (2011); Lin et PID al.designs, (1990); Bing etiscontrol al. (2015)). A to major drawback of such however, their inability A particularly interesting technology hyperthermia is hybrid sound (MR-HIFU), which enables heatfor delivery to internal magnetic-resonance-guided high-intensity focused ultradrawback of such designs, however, isal. their inability to take sound (MR-HIFU), which enables heat delivery to internal (2011); Lin et al. (1990); Bing et (2015)). A major the body’s future thermal behavior and the restrictive acdrawback of such designs, however, is their inability to take magnetic-resonance-guided high-intensity focused ultratissues(MR-HIFU), with spatial accuracy in theheat millimeter range, com- the body’s future thermal behavior and the restrictive acsound which enables delivery to internal tissues with spatial in millimeter range, comoffuture such designs, however, isand their inability take tuator constraints into account when computing theto input, the body’s thermal behavior the restrictive acsound (MR-HIFU), which enables heat toresponse internal tissueswith with real-time spatial accuracy accuracy in the theof millimeter range, com- drawback bined monitoring thedelivery thermal tuator constraints into account whenand computing the input, bined with real-time monitoring of the thermal response body’s future affects thermal behavior theand restrictive acwhich negatively treatment duration. tuator constraints into account whenquality computing the input, tissues with real-time spatial accuracy in theofmillimeter range, com- the bined with monitoring the thermal response which treatment and ⋆ This research tuator constraintsaffects into account whenquality computing the input, which negatively negatively affects treatment quality and duration. duration. has been made possible of by the Cancer Society bined with real-time monitoring theDutch thermal response ⋆ This research has been made possible by the Dutch Cancer Society In thisnegatively respect, model predictive quality control (MPC) offers ⋆ This which affects treatment duration. and theresearch Netherlands Organisation for by Scientific Research has been made possible the Dutch Cancer (NWO) Society In this respect, model predictive control and (MPC) offers and the Netherlands Organisation for Scientific Research (NWO) superior potential for hyperthermia. This has been recIn this respect, model predictive control (MPC) offers ⋆ as This part of theirhasjoint Partnership Programme: “Technology for been made possible by the Dutch Cancer (NWO) Society and theresearch Netherlands Organisation for Scientific Research superior potential for hyperthermia. This has been recas part of their joint Partnership Programme: “Technology for In this by respect, model predictive control (MPC) superior potential for hyperthermia. Thisresults has been recognized some researchers, of which the areoffers pubOncology”. This is partially financed by the PPP Allowance as part theirproject jointOrganisation Partnership Programme: “Technology for and the of Netherlands for Scientific Research (NWO) ognized by some researchers, of which the results are pubOncology”. This project is partially financed by the PPP Allowance superior potential for hyperthermia. This has been reclished in, for example, Arora et al. (2007) for the control ognized by some researchers, of which the results are pubmade available by Top Sector Life Sciences & Health. This work is Oncology”. is partially financed by the“Technology PPP Allowance as part of This theirproject joint Partnership Programme: for lished in, for example, Arora et al. (2007) for the control of of made available by Top Sector Life Sciences & Health. This work is ognized researchers, of al. which thefor results are pubalso partially funded by Sector the European Union & via IPaCT Project. lished in,byforsome example, Arora et (2007) the control of made available byproject Top Life Sciences This work is Oncology”. This is partially financed byHealth. the PPP Allowance also partially funded by the European Union via the IPaCT Project. lished in, for example, Arora et al. (2007) for the control of also partially funded by Sector the European Union & viaHealth. the IPaCT made available by Top Life Sciences This Project. work is also partially funded by the European Union via the IPaCT Project. 2405-8963 © 2018 2018, IFAC (International Federation of Automatic Control) Copyright © IFAC 223 Hosting by Elsevier Ltd. All rights reserved. Copyright © under 2018 IFAC 223 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 223 10.1016/j.ifacol.2018.11.012 Copyright © 2018 IFAC 223
2018 IFAC NMPC 192 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
thermal dose in a single point based on a one-dimensional temperature distribution, De Bever et al. (2014) for thermal dose control in two-dimensional systems using heavily simplified models and a fixed sonication trajectory, and more recently in Sebeke et al. (2017) for temperature control in two-dimensional systems with online adaption of the heating location and power. Although the initial results of online optimization-based control in the aforementioned works look promising, it is also clear that the full potential of MPC is far from reached. More sophisticated MPC setups specifically tailored to hyperthermia can fulfill the envisioned potential, resulting in more reliable, reproducible and highperforming 3D temperature control, which in turn may serve as an enabler for widespread and successful clinical use of MR-HIFU for hyperthermia. It is the objective of this paper to provide a significant step forward in the development of this societally highly relevant technology. To this end, an offset-free MPC design is developed for planar heating using an MR-HIFU system. A Luenbergertype observer is employed to improve the MR thermometry measurements and to construct a disturbance estimate. The latter enables the controller to compensate for the effects of constant (or slowly varying) model errors, allowing for enhanced control performance in the clinic, where such model discrepancies are inevitable. Simulations are performed to demonstrate the developed model-based algorithm’s potential in terms of achieving desired steadystate behavior despite severe plant-model mismatch. The remainder of this paper is organized as follows. In Section 2 relevant aspects related to the clinical application of mild hyperthermia are discussed in more detail to clarify the challenges involved with feedback control for this application, and to motivate the desire to develop MPC for MR-HIFU. Next, Section 3 explains how the thermal model is obtained and Section 4 covers the controller design procedure. Simulations are performed to evaluate the algorithm’s ability of removing the steady-state offset otherwise induced by model mismatch, the results of which are analyzed in Section 5. Finally, Section 6 summarizes the key achievements and corresponding observations. 2. PROBLEM DESCRIPTION This section introduces some practical aspects concerning the hyperthermia treatment and setup, and discusses the motivation for using MPC in such an application. 2.1 MR-HIFU hyperthermia setup Although the MPC setup we propose is generic in nature, in this paper it is designed for the Philips Sonalleve shown in Fig. 1, which is a commercial MR-HIFU system already being used in clinics to treat uterine fibroids and for bone pain palliation therapy. It is a combined setup of a standard MRI scanner used for noninvasive thermography and a dedicated trolley-tabletop in which an MR-compatible HIFU transducer is integrated. 2.2 HIFU applicator This system uses a phased-array HIFU transducer to generate the ultrasound waves. It consists of 256 elements 224
Fig. 1. Philips Sonalleve MR-HIFU therapy platform. focal plane focal point
transducer driving signals
tumor
patient table transducer axis
Fig. 2. Schematic of HIFU beam into focal plane in tumor area, with focal point by electronic beam steering. of which the phases and amplitudes can individually be chosen such that by interference a focal spot is created. By modulating these settings, which is referred to as electronic beam steering, the lateral location of the focal spot can be chosen anywhere on the focal plane within an 18 mm diameter circle around the transducer axis, as depicted in Fig. 2. Since the size of the focal spot is typically smaller than the tumor, it must be scanned through the treatment volume using beam steering or by moving the transducer. In this paper, we focus on beam steering in the focal plane using a predetermined set of discrete possible focus locations. 2.3 Motivation for offset-free MPC in hyperthermia To fully benefit from the hyperthermia-induced effects such as increased blood flow and oxygenation, it is crucial to maintain temperatures above 41 ◦ C in the region of interest (ROI), typically the tumor and possibly some adjacent tissue, during the entire 90 minute treatment. On the other hand, temperatures above 43 ◦ C are undesired since they may cause reverse effects such as vascular shutdown (Hijnen et al. (2014)). Furthermore, temperature elevation outside the ROI must be avoided to prevent damage to healthy tissue. MPC is considered superior to PID-based strategies for application in MR-HIFU hyperthermia, as it is able to reduce treatment time by exploiting beneficial (future) behavior, e.g., heat diffusion via conduction, and can explicitly incorporate actuator constraints such as the inability to actively remove heat from inside the body using HIFU. Unfortunately, model-based strategies are inherently accompanied by the possibility of modeling errors, which in the thermal models for hyperthermia treatment are highly likely to occur due to the high variability of spatially distributed and time/temperaturevarying tissue parameters. Inadequately accounting for such mismatches may result in insufficient heating of the target region or in the overheating of healthy tissue, thereby significantly deteriorating treatment quality.
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
acoustic properties of the intermediate tissue. In this paper, however, F is modeled by a radially symmetric twodimensional Gaussian distribution centered around the focus location rf with standard deviation σf = 2.4 mm �r − rf (t)�2 α , (3) exp − F (r, t) = 2πσf2 2σf2
MR-HIFU setup and patient
Desired acoustic power
with r ∈ Ω, t ∈ R≥0 , and scaling factor α ∈ R>0 .
MRT data MPC
193
3.2 Discrete-time state-space model
Model-based observer T setpoint, constraints
Fig. 3. Block scheme of envisioned fully automated MRHIFU feedback loop for patient treatment. Given this situation, we propose an MPC setup novel for hyperthermia inspired by an offset-free MPC algorithm as in Pannocchia and Rawlings (2003). Fig. 3 depicts the resulting feedback scheme, where the observer constructs the model-based temperature and disturbance estimates from the noninvasively obtained MR thermometry data, such that the MPC scheme can compute the optimal power distribution to be generated by the HIFU actuator. 3. THERMAL PROCESS MODELING In this section, a model of the body’s thermal response will be discussed, from which the state-space models used by the observer and MPC algorithm will be derived. 3.1 Bioheat model The tissue’s temperature change in response to the applied heating is modeled using Pennes’ bioheat equation, see Pennes (1948), given by ∂T (r, t) = k∇2 T (r, t) + Q(r, t) − wb (r)cb (T (r, t) − Tb ), ρc ∂t (1) where T : Ω × R≥0 → R is the temperature profile, such that T (r, t) denotes the temperature at time t ∈ R≥0 and location r = [rx , ry ]⊤ ∈ Ω ⊂ R2 with Ω being the patient domain in the focal plane. The parameters ρ, c, k ∈ R>0 denote the tissue’s volumetric mass density, specific heat capacity, and thermal conductivity, and wb : Ω → R≥0 and cb ∈ R>0 the blood’s perfusion coefficient and specific heat capacity, respectively. In this paper, the tissue is assumed to be homogeneous with respect to most parameters (ρ, k, c, and cb ), but not the blood perfusion rate wb (r), which is assumed significant only in tissue in the direct vicinity of blood vessels, and zero elsewhere. Note, however, that our MPC setup can be extended to incorporate the fully inhomogeneous case as well. The power deposition density described by Q : Ω × R≥0 → R≥0 depends on the acoustic deposition intensity profile F : Ω × R≥0 → R≥0 and linearly scales with the sonication power P : R≥0 → R≥0 according to Q(r, t) = F (r, t)P (t), r ∈ Ω, t ∈ R≥0 . (2) In reality, for a given focus location rf (t) at some time t, the intensity F (r, t) also depends on factors such as the 225
Assuming zero out-of-plane interaction, (1)-(3) are discretized spatially using the central difference scheme on a 39 × 39 grid with a voxel size of 2 × 2 mm2 and temporally using forward Euler with a sampling interval of Ts = 2.5 s, resulting in the discrete-time state-space model xk+1 = Ap xk + Bp uk , (4a) (4b) yk = xk + vk , with k ∈ N connecting to real time tk = kTs and the states xk ∈ Rn , n = 392 = 1521 representing the voxel temperatures. We select 48 voxels within the ROI, at the center of which we allow sonication by steering the focal point rf to this location. These locations are referred to as sonication points. The input uk ∈ Rm with m = 48 is the applied sonication power averaged over the sampling interval at each of these sonication points. Each individual column of Bp ∈ Rn×m captures the spatially discretized power deposition as in (2) and (3) corresponding to rf at a single sonication point, and its temporally discretized relation to voxel temperature change according to (1). The spatial discretization is chosen in such a manner that the centers of the voxels coincide with the points measured by MR thermometry. Also, the sampling interval Ts is larger than the time required for measurement, implying that at each discrete time k ∈ N a full new temperature map is available and thus the output yk ∈ Rn can indeed be modeled as a full (noisy) state measurement as in (4b). In practice, the MR data is corrupted by noise, which for the purpose of this paper can be approximated by uncorrelated zero-mean Gaussian white vk ∈ Rn , √ noise ◦ k ∈ N, with standard deviation σv = 0.2 C. 4. CONTROLLER DESIGN In this section, the state-space models used by the controller and observer are given first. Then, the control objectives as described in Subsection 2.3 are formulated as a constrained optimization problem. 4.1 Prediction model To enable the controller to compensate for the steadystate offset induced by model mismatch or slowly varying disturbances, integral action is incorporated using disturbance estimation as described in Pannocchia and Rawlings (2003). The state-space representation of the prediction model is given by xi+1|k = Axi|k + Bui|k + di|k , (5a) di+1|k = di|k , (5b) yi|k = xi|k , (5c) where xi|k , ui|k , yi|k ∈ Rn denote the predicted states, inputs, and outputs, respectively, at i ∈ N time steps ahead
2018 IFAC NMPC 194 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
of the prediction sequence’s starting time k. Furthermore, di|k ∈ Rn denotes the virtual state disturbance, which is assumed to be constant over the prediction horizon. Note that we allow model mismatch, i.e., A �= Ap and B �= Bp , and that we did not incorporate measurement noise in the prediction model since E(vk ) = 0. As proposed in Theorem 1 of Pannocchia and Rawlings (2003), we first verify that the nominal model (i.e., the model without additional disturbance di|k ) is detectable (in fact, observable). Second, we evaluate the rank condition stated in the theorem, which using our disturbance model becomes I − A −I rank = 2n, (6) I 0
to conclude that (5) is detectable, and therefore an asymptotically stable estimator exists given the proposed disturbance model. This property also guarantees that for a stable closed-loop system, a feasible steady-state setpoint, and constant disturbances, there will be zero offset between the realized and desired temperature distributions in case no constraints are active at steady state. 4.2 State and disturbance estimation To facilitate this offset-free implementation, and to simultaneously obtain improved temperature estimates with respect to the noise-corrupted MR thermometry readings, we design the state and disturbance estimator x ˆk = Aˆ xk−1 + Buk−1 + dˆk−1 + Lx (yk − yˆ− ), (7a) k
(7b) dˆk = dˆk−1 + Ld (yk − yˆk− ), where yˆk− = Aˆ xk−1 + Buk−1 + dˆk−1 (8) denotes the model-based measurement estimate at time k before applying the correction by output injection to obtain x ˆk . Note that we choose the observer model (7) equal to the prediction model (5). Since we have full (noisy) state measurement, the system (5) is observable, allowing for arbitrary pole placement in (7) via Lx and Ld . We choose the diagonal matrices 1 1 In , Lx = In , Ld = (9) 3 30 which cause the state and disturbance estimation errors to converge to within 5% of their initial value in approximately 22.5 and 225 seconds, respectively. 4.3 Optimization problem
T [◦ C] 43 42 41 40
T +ǫ T −ǫ
Transducer axis ǫ
ǫ
37
R
r2Ω
18 mm
S 32 mm
Fig. 4. Schematic cross-section of the temperature objectives for the circular regions R and S centered at the transducer axis. The maximum violations ǫ and ǫ are shown for some over- and underheated temperature distributions such that T1 ≤ T + ǫ and T2 ≥ T − ǫ. lower value on Ω \ S to safeguard healthy tissue. Within the aforementioned temperature range, preference is given to a flat temperature distribution in the ROI, which is described by the reference temperature Tr . To formulate these objectives mathematically, let us denote the elements of x corresponding to a point within R by xR = Hx ∈ RnR , where H ∈ {0, 1}nR×n is a matrix with exactly one 1 in each row. Furthermore, let xr ∈ RnR and x, x ∈ Rn denote the voxel-wise temperature reference and lower/upper bounds corresponding to the values of Tr , T , and T on the discrete voxel locations, respectively. The voxels’ temperature deviations with respect to the reference in the ROI are then given by xR − xr . The maximum violations of the temperature bounds follow from ǫ = �max(x − x, 0)�∞ ≥ 0 and ǫ = �max(x − x, 0)�∞ ≥ 0, where the maximum operator is used component-wise, and are collected in ǫ = [ǫ, ǫ]⊤ ∈ R2≥0 . The control objective is now formulated as the constrained optimization problem N −1 min ℓ(xi|k , ui|k , ǫi|k ) + ℓ(xN |k , 0, ǫN |k ) (10a) u , i|k
i∈[0,N −1]
i=0
with, for some i, k ∈ N, the stage cost r r ⊤ R ℓ(xi|k , ui|k , ǫi|k ) = (xR i|k − x ) Q(xi|k − x ) ⊤ + u⊤ i|k Rui|k + fǫ ǫi|k ,
The temperature objectives have been discussed in Subsection 2.3 and are schematically depicted in Fig. 4 in cross-section perspective. T : Ω → R and T : Ω → R represent the (location-dependent) lower and upper temperature bounds, respectively, defining the desired temperature range (green) over the entire patient domain Ω. Underand overheating occurs for temperatures within the (light) blue and (light) red areas in the figure, respectively. The concentric circular regions R and S denote the ROI and the region additionally including the tissue surrounding the ROI, i.e., R ⊂ S ⊂ Ω. T is nonzero only within R where sufficient heating is desired . T features an elevated plateau on S, preventing reversal of the temperaturedependent mechanisms due to overheating, and has a 226
Tr T;T T1 T2
subject to x − 1ǫi|k xi+1|k di+1|k ui|k
(10b)
i ∈ [0, N ], (10c) ≤ xi|k ≤ x + 1ǫi|k , = Axi|k + Bui|k + di|k , i ∈ [0, N ], (10d) = di|k , i ∈ [0, N − 1], (10e) ≥ 0, i ∈ [0, N − 1], (10f)
1⊤ ui|k ≤ utot ,
i ∈ [0, N − 1]. (10g)
Inequality (10c), in which 1 denotes an all-ones vector of suitable dimension, captures the temperature bounds as soft constraints (note that ǫi|k = [ǫi|k , ǫi|k ]⊤ ), (10d)(10e) ensure satisfaction of the dynamics according to (5) with the prediction sequence’s initial conditions given by ˆk and d0|k = dˆk , and finally (10f)-(10g) describe x0|k = x
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
which are normalized with respect to the number of corresponding variables for more intuitive balancing of the objectives’ relative contribution to the cost function. The critical importance of achieving temperatures within the desired range is expressed by the linear term (which is possible since ǫ ≥ 0 by definition) with high weighting. Using only the maximum violation, as opposed to penalizing each voxel individually, improves temperature homogeneity during the initial heat-up phase. To also encourage a flat temperature distribution (within the temperature bounds) in the ROI at steady steady, a quadratic term with relatively lower cost is added to penalize this tracking error. In addition to the aforementioned control objectives, a small weight on the input power is incorporated to favor solutions that involve less heating.
195
43 Temp. in ROI [◦ C]
the actuator constraints with utot = 60 W being the maximum allowable power being applied to the treated region. Furthermore, a horizon of N = 5 is used, and the weights are chosen to be 0.01 10 100 I, fǫ = I, R = , (11) Q= 200 nR m
42 41 xR avg R xR min , xmax R x ˆavg x ˆR ˆR max min , x
40 39 38 37
0
100
200
300 Time [s]
400
500
600
400
500
600
(a)
43 Temp. in ROI [◦ C]
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
42 41 40 39 38 37
0
100
200
5. PERFORMANCE ANALYSIS
300 Time [s] (b)
Next, we perform a simulation in which using the nominal model results in temperature overestimation. To this end, the thermal conductivity k in (1) is intentionally underestimated and the power deposition intensity’s scaling factor α in (3) is overestimated in the prediction and estimation models, both by a factor of 5, affecting A and B, respectively, in (5) and (7). Additionally, the observer and controller models are set to assume a homogeneous perfusion rate wb (r) = wb , thereby neglecting the local influence of blood vessels present in the plant. The results are shown in Fig. 5b. Initially, the disturbance estimator is disabled, i.e., dˆk = 0 and Ld = 0 for k ∈ N corresponding to tk < 300 s. This causes a mismatch between x ˆR and xR , resulting in a downward offset in the realized xR with respect to its reference xr . From tk = 300 s (vertical dashed black line) onwards, the offset-free algorithm is active by setting Ld as ˆ the achieved in (9). At tk = 600 s, after convergence of d, temperature distribution coincides with that of Fig. 5a, illustrating how offset-free MPC can significantly improve hyperthermia treatment quality under model error. This is more clearly illustrated in Fig. 6, where similarly to Fig. 4 cross-sections of the temperature objectives and distributions through the center of R (at rx = 0.128 m) 227
Fig. 5. The average (solid) and lowest/highest (dashed) voxel temperatures in the ROI of the plant xR (black) and observer x ˆR (gray) in case of (a) no plant-model mismatch, and (b) significant model error where disturbance estimation is enabled at t = 300 s. Temperature [◦ C]
A simulation study was carried out to verify the effectiveness of the proposed control algorithm in terms of compensating for the effects of plant-model mismatch for the MR-HIFU system under consideration. For reference, a simulation without model error is performed first, i.e., where A = Ap and B = Bp , of which the results are shown in Fig. 5a. Initially, the entire patient volume is assumed to be 37 ◦ C. The estimated temperature x ˆR (gray) closely R resembles its true counterpart x (black) in the ROI, as is shown in the figure by their lowest/highest (dashed) and average value (solid), enabling the controller to accurately heat the target region in a satisfactory manner. The current MPC’s worst-case computation time is found to be 1 s, indicating it is sufficiently fast given Ts = 2.5 s.
43 42 41 40 39 38 37 0.09 0.10
(a) (b) (c) Fi
0.11
0.12
0.13 ry [m]
0.14
0.15
0.16
Fig. 6. Voxel temperature xk at cross-section through center of ROI without plant-model mismatch at tk = 600 s (a), and with mismatch at tk = 300 s (b) and tk = 600 s (c), alongside the acoustic intensity around in-plane sonication points Fi . The gray areas indicate the cross-section voxels with increased perfusion. are shown, alongside the (scaled) acoustic intensities Fi = F (rf,i ), i = 1, . . . , 6, from (3) with center locations rf,i at the six in-plane sonication points. The temperature profile corresponding to Fig. 5a at tk = 600 s (black) shows the MPC’s optimal performance without model error. It is within the desired temperature range, i.e., x ≤ x ≤ x, and tends to xr in R, however due to the system restrictions it is not a perfectly flat distribution. That is, the acoustic intensity profiles (blue) do not admit further heating of the areas with locally increased perfusion (gray) without causing additional overheating outside these areas. As previously found in Fig. 5b at tk = 300 s, Fig. 6 shows that while introducing plantmodel mismatch with disturbance estimation disabled
ry [m]
2018 IFAC NMPC 196 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
0.16
0.6
0.15
0.4
0.14
0.2
0.13
0
0.12
−0.2
0.11 0.10
−0.4
0.09
−0.6 0.10
0.14 0.12 rx [m]
0.16
Fig. 7. Disturbance estimate dˆ at time tk = 600 s, with additionally the edges of the circular R (solid black) and S (dashed black), sonication points (black dots), and voxels with increased perfusion (red squares). causes insufficient heating (solid red), optimal performance is recovered by the offset-free algorithm after convergence at tk = 600 s (dashed red). The disturbance estimate after convergence at tk = 600 s is visualized in Fig. 7. Indeed, overestimation occurs in the entire ROI due to the mismatch in k and α, which ˆ Moreover, the is compensated for by negative values in d. figure shows that the disturbance estimator correctly identifies the voxels in which the locally increased perfusion is neglected (red outline), and thus is able to counteract the corresponding effects. In the figure, some voxels outside S have a slightly positive disturbance estimate, which is a result of the propagation of measurement noise into dˆ via the output injection in (7). 6. CONCLUSION In this paper, an offset-free MPC scheme is presented to enhance control performance of an MR-HIFU hyperthermia system used in cancer treatment. The proposed design achieves more accurate 2D temperature control than the techniques currently used in clinics. Specifically, the model-based algorithm is able to explicitly take into account (a) the restrictive input constraints and (b) the body’s future thermal response, thereby enabling more optimal heating and shorter treatment times. Using disturbance estimation, the control design eliminates the steady-state offset otherwise resulting from modeling error, thus enhancing hyperthermia treatment quality. The MPC setup’s effectiveness in terms of compensating for unmodeled perfusion and inaccurate conduction/heating parameters is demonstrated using simulation studies. REFERENCES Arora, D., Skliar, M., Cooley, D., and Roemer, R.B. (2007). Constrained Predictive Control of Thermal Therapies for Minimum-Time Delivery of Thermal Dose. IEEE Trans. Control Syst. Technol., 15(6), 1030– 1037. Bing, C., Nofiele, J., Staruch, R., et al. (2015). Localised hyperthermia in rodent models using an MRIcompatible high-intensity focused ultrasound system. Int. J. Hyperth., 31(8), 813–822. 228
Datta, N., Ord´on ˜ ez, S.G., Gaipl, U., et al. (2015). Local hyperthermia combined with radiotherapy and-/or chemotherapy: Recent advances and promises for the future. Cancer Treat. Rev., 41(9), 742–753. De Bever, J., Todd, N., Payne, A., Christensen, D.A., and Roemer, R.B. (2014). Adaptive model-predictive controller for magnetic resonance guided focused ultrasound therapy. Int. J. Hyperth., 30(7), 456–470. Enholm, J., Kohler, M., Quesson, B., et al. (2010). Improved Volumetric MR-HIFU Ablation by Robust Binary Feedback Control. IEEE Trans. Biomed. Eng., 57(1), 103–113. Franckena, M., Fatehi, D., De Bruijne, M., et al. (2009). Hyperthermia dose-effect relationship in 420 patients with cervical cancer treated with combined radiotherapy and hyperthermia. Eur. J. Cancer, 45(11), 1969–1978. Hijnen, N., Langereis, S., and Gr¨ ull, H. (2014). Magnetic resonance guided high-intensity focused ultrasound for image-guided temperature-induced drug delivery. Adv. Drug Deliv. Rev., 72, 65–81. Issels, R.D., Lindner, L.H., Verweij, J., et al. (2018). Effect of Neoadjuvant Chemotherapy Plus Regional Hyperthermia on Long-term Outcomes Among Patients With Localized High-Risk Soft Tissue Sarcoma. JAMA Oncol., 4(4), 483–492. Lin, W.L., Roemer, R.B., and Hynynen, K. (1990). Theoretical and experimental evaluation of a temperature controller for scanned focused ultrasound hyperthermia. Med. Phys., 17(4), 615–625. Mallory, M., Gogineni, E., Jones, G.C., Greer, L., and Simone, C.B. (2016). Therapeutic hyperthermia: The old, the new, and the upcoming. Crit. Rev. Oncol. Hematol., 97, 56–64. Maloney, E. and Hwang, J.H. (2015). Emerging HIFU applications in cancer therapy. Int. J. Hyperth., 31(3), 302–309. Mougenot, C., Quesson, B., De Senneville, B.D., et al. (2009). Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU). Magn. Reson. Med., 61(3), 603–614. Pannocchia, G. and Rawlings, J.B. (2003). Disturbance models for offset-free model-predictive control. AIChE J., 49(2), 426–437. Partanen, A., Yarmolenko, P.S., Viitala, A., et al. (2012). Mild hyperthermia with magnetic resonance-guided high-intensity focused ultrasound for applications in drug delivery. Int. J. Hyperth., 28(4), 320–336. Pennes, H. (1948). Analysis of Tissue and Arterial Blood Temperatures in the Resting Human Forearm. J. Appl. Physiol., 1(2), 93–122. Sebeke, L., Luo, X., Heijman, E., et al. (2017). Predictionbased controller for MR-HIFU mediated hyperthermia. In Eur. Symp. Focus. Ultrasound Ther., 23. Staruch, R., Chopra, R., and Hynynen, K. (2011). Localised drug release using MRI-controlled focused ultrasound hyperthermia. Int. J. Hyperth., 27(2), 156–171. Staruch, R.M., Hynynen, K., and Chopra, R. (2015). Hyperthermia-mediated doxorubicin release from thermosensitive liposomes using MR-HIFU: Therapeutic effect in rabbit Vx2 tumours. Int. J. Hyperth., 31(2), 118–133.
ScienceDirect
IFAC PapersOnLine 51-20 (2018) 191–196
Offset-free model predictive control for Offset-free model predictive control for Offset-free model predictive control for enhancing MR-HIFU hyperthermia in Offset-free model predictive control for enhancing MR-HIFU hyperthermia in enhancingcancer MR-HIFU hyperthermia in ⋆ treatment ⋆⋆ enhancingcancer MR-HIFU hyperthermia in treatment cancer treatment ⋆ cancer ∗ ∗∗ ∗ D.A. Deenen ∗∗ E. Maljaarstreatment ∗ L. Sebeke ∗∗ Bram de Jager ∗
D.A. Deenen ∗ E.∗∗∗ Maljaars ∗ L. ∗∗ Sebeke ∗∗ Bram de Jager ∗ ∗ Heijman H. Gr¨ ullL. W.P.M.H.Bram Heemels D.A. E. Deenen E.∗∗∗ Maljaars de Jager ∗∗ Sebeke ∗ E. Heijman H. Gr¨ u Heemels ∗ ∗ll ∗∗ W.P.M.H. ∗∗ ∗ ∗∗∗ ∗ D.A. E. Deenen E. Maljaars Bram de Jager Heijman H. Gr¨ ullL. Sebeke W.P.M.H. Heemels ∗∗∗ ∗∗ ∗ ∗ E. Heijman H. Gr¨ ull W.P.M.H. Heemels ∗ Eindhoven University of Technology, Mechanical Engineering, University of Technology, Mechanical Engineering, ∗ Eindhoven Control Systems Technology, Eindhoven, The Netherlands Eindhoven University of Technology, Mechanical Engineering, ∗ Control Systems Technology, Eindhoven, The Netherlands Eindhoven University of [email protected]) Technology, Mechanical Engineering, (e-mail: Control Systems Technology, Eindhoven, The Netherlands (e-mail: [email protected]) ∗∗ Control Systems Technology, Eindhoven, TheGermany Netherlands (e-mail: [email protected]) Hospital Cologne, Cologne, ∗∗ University University Hospital Cologne, Cologne, Germany ∗∗ ∗∗∗ [email protected]) Philips (e-mail: Research Germany, University Hospital Cologne, Aachen, Cologne, Germany Germany ∗∗∗ Philips Research Germany, ∗∗ ∗∗∗ University HospitalGermany, Cologne, Aachen, Cologne, Germany Germany Philips Research Aachen, Germany ∗∗∗ Philips Research Germany, Aachen, Germany Abstract: Abstract: In In this this paper, paper, an an offset-free offset-free MPC MPC scheme scheme is is developed developed to to improve improve performance performance and robustness of temperature control inMPC magnetic focused Abstract: In this paper, an offset-free schemeresonance-guided is developed to high-intensity improve performance and robustness of temperature control in magnetic resonance-guided high-intensity focused Abstract: In this paper,cancer an offset-free is design developed to high-intensity improve and robustness of temperature control inMPC magnetic resonance-guided focused ultrasound hyperthermia treatments. The scheme proposed effectively handlesperformance constraints ultrasound hyperthermia cancer treatments. The proposed design effectively handles constraints robustness of temperature control to in achieve magnetic resonance-guided focused and exploits the body’s thermal faster and moreeffectively evenlyhigh-intensity distributed heating of ultrasound hyperthermia cancerbehavior treatments. The proposed design handles constraints heating of and exploitshyperthermia the body’s thermal behavior to achieve faster and moreeffectively evenly distributed ultrasound cancertreatment treatments. The proposed design handles constraints the tumor leading to improved with respect to the methods currently applied and exploits the body’s thermal behavior toquality achieve faster and more evenly distributed heating of the tumor leading to improved treatment with respect to the methods currently applied and exploits the proposed body’s thermal behavior toquality achieve faster andestimation more evenly heating of in The scheme incorporates disturbance todistributed remove the steadytheclinics. tumor leading to improved treatment quality with respect to the methods currently applied in clinics. The proposed scheme incorporates disturbance estimation to remove the steadythe tumor to from improved treatment qualitydisturbance with respect to given the methods currently in clinics. The proposed scheme incorporates estimation to remove the applied steadystate offsetleading resulting plant-model mismatch, which for the healthcare application is state offsetThe resulting from scheme plant-model mismatch, which for the given healthcare is in clinics. proposed incorporates disturbance estimation to removeapplication the steadyinevitable practice due the diversity and Simulation state offsetin mismatch, which for in thepatients given healthcare application is inevitable inresulting practice from due to toplant-model the large large physiological physiological diversity in patients and tumors. tumors. Simulation state offset resulting from mismatch, which for in thepatients healthcare application is results are in presented demonstrate algorithm’s capability togiven improve performance. inevitable practice to due toplant-model the largethis physiological diversity androbust tumors. Simulation results are presented demonstrate algorithm’s capability to improve performance. inevitable practice to due to the largethis physiological diversity in patients androbust tumors. Simulation results are in presented to demonstrate this algorithm’s capability to improve robust performance. © 2018,are IFAC (International Federation ofthis Automatic Control) Hosting to by improve Elsevier Ltd. All rights reserved. results presented to demonstrate algorithm’s capability robust performance. Keywords: Keywords: Healthcare, Healthcare, hyperthermia, hyperthermia, offset-free offset-free MPC, MPC, cancer cancer treatment, treatment, high-intensity high-intensity Keywords: Healthcare, hyperthermia, offset-free MPC, cancer treatment, high-intensity focused ultrasound. focused ultrasound. Keywords: Healthcare, hyperthermia, offset-free MPC, cancer treatment, high-intensity focused ultrasound. focused 1. ultrasound. INTRODUCTION via volumetric temperature maps obtained using an MRI 1. via volumetric temperature maps using an 1. INTRODUCTION INTRODUCTION via volumetric temperature maps obtained obtained usingallows an MRI MRI scanner, see Maloney and Hwang (2015). This for scanner, see Maloney and Hwang (2015). This allows for 1. INTRODUCTION via volumetric temperature maps obtained using an MRI a completely noninvasive procedure, thereby significantly scanner, see Maloney and Hwang (2015). This allows for Mild local hyperthermia involves the heating of a specific a completely noninvasive procedure, thereby significantly Mild local hyperthermia involves the heating of a specific see Maloney Hwang (2015). This allows the for to patientand well-being. However, realizing acontributing completely noninvasive procedure, thereby significantly Mild local hyperthermia the heating of a specific target volume inside the involves human body to temperatures of scanner, contributing to patient well-being. However, the target inside human body to temperatures of ◦volume completely noninvasive procedure, thereby realizing significantly desired temperature distribution and accurately maintainto patient well-being. However, realizing the Mild local hyperthermia the heating of a specific target inside the the involves human body toNumerous temperatures of acontributing C for approximately 90 minutes. clinical 39-45 ◦volume temperature distribution and accurately maintain39-45 ◦ C for approximately 90 minutes. Numerous clinical desired to course patient However, realizing the desired temperature distribution and accurately ing it over the ofwell-being. an entire treatment in maintaina clinical targethave volume inside human bodyisto temperatures of contributing 39-45 trials thatthe hyperthermia aNumerous potent adjuvant C forshown approximately 90 minutes. clinical ing it over the course of an entire treatment in a clinical trials have shown that hyperthermia is a potent adjuvant ◦ desired temperature distribution and accurately maintainsetting is no simple task. Combined with the fact that ing it over the course of an entire treatment in a clinical 39-45 C for approximately 90 minutes. Numerous clinical therapy for shown cancerthat treatment, see for example Datta et al. setting is no simple task. Combined with the fact that trials have hyperthermia is a potent adjuvant therapy for shown cancerthat treatment, see for example Datta et al. ing over thesimple coursetask. entire therapies, treatment in fact a of clinical the itbeneficial effects of an thermal and mild is no Combined with the that trials hyperthermia isal.a (2018). potent adjuvant (2015);have Mallory et al. (2016); Issels It signiftherapy for cancer treatment, see foretexample Datta et al. setting the beneficial effects task. of thermal therapies, andfact of mild (2015); Mallory et al. (2016); Issels et al. (2018). It signifsetting is no in simple Combined with the that hyperthermia particular, are closely related to the the beneficial effects of thermal therapies, and oftissue mild therapy for cancer treatment, seeprimary foretexample Datta et al. icantly increases the efficacy of methods such as (2015); Mallory et al. (2016); Issels al. (2018). It signifin particular, are closely related to the icantly increases the efficacy of primary methods such as hyperthermia beneficialtruly of thermal therapies, oftissue mild hyperthermia ineffects particular, are closely related and tosee the tissue temperatures achieved during treatment, Franck(2015); Mallory etthe al.efficacy (2016); Issels et al. (2018). signifchemo- increases and radiotherapy by sensitizing the treatedItsuch tumoricantly of primary methods as the temperatures truly achievedare during treatment, Franckchemoand radiotherapy by sensitizing the treated tumorin particular, closely related tosee the tissue ena this strongly motivates the temperatures truly achieved during treatment, see Franckicantly increases efficacy of methods as hyperthermia chemoandAlso, radiotherapy by the treated such tumorous tissue. itthe enables thesensitizing useprimary of temperature-sensitive ena et et al. al. (2009), (2009), this strongly motivates the development development ous tissue. Also, it enables the use of temperature-sensitive trulythis achieved during treatment, see Franckof feedback controllers for MR-HIFU. ena et al. (2009), strongly motivates the development chemoandAlso, radiotherapy the treated tumorliposomes for heat-mediated drug release, allowing for temperatures ous tissue. it enablesby thesensitizing use of temperature-sensitive feedback controllers for MR-HIFU. liposomes for heat-mediated drug release, allowing for of et al. (2009), this strongly motivates the development of feedback controllers for MR-HIFU. ous Also, it enables the anticancer use of temperature-sensitive moretissue. efficient of the drugsallowing specifically liposomes for delivery heat-mediated drug release, for ena Currently, (experimental) clinical more efficient delivery of the anticancer drugs specifically feedback (experimental) controllers for MR-HIFU. Currently, clinical applications applications of of local local hyhyliposomes forregion heat-mediated drug release, allowing for of to theefficient tumor at lower drug concentrations, more delivery of thesystemic anticancer drugs specifically perthermia(experimental) by MR-HIFU typically involve a (mostly) Currently, clinical applications of localprehyto the tumor region at lower systemic drug concentrations, perthermia by MR-HIFU typically involve a (mostly) premore efficient delivery of the anticancer drugs specifically to the tumor region at lower systemic drug concentrations, as in Hijnen et al. (2014); Staruch et al. (2015). Most Currently, clinical applications of local hydetermined(experimental) sonication plan, possibly extended with simple perthermia by MR-HIFU typically involve a (mostly) preas in Hijnen et al. (2014); Staruch et al. (2015). Most sonication plan, possiblyinvolve extended with simple to the tumor the region lower drug importantly, showsystemic that local does determined as in Hijnen et results al. at(2014); Staruch ethyperthermia al.concentrations, (2015). Most perthermia by MR-HIFU typically a (mostly) preonline MR-driven feedback algorithms. Examples include determined sonication plan, possibly extended with simple importantly, show Staruch that localethyperthermia does online MR-driven feedback algorithms. Examples include as Hijnenthe et results al. (2014); al. (2015). Most notinincrease toxicity inshow healthy tissuehyperthermia and the thereby importantly, the results that local does determined sonication plan, algorithms. possibly extended with and/or simple binary MR-driven strategies scale power online feedback Examples not increase toxicity in healthy tissue and the thereby strategies to to scale the the sonication sonication power include and/or importantly, results that local hyperthermia does induced side the effects, making it highly clin- binary not increase toxicity inshow healthy tissueappealing and the for thereby online feedback algorithms. Examples include length MR-driven of the heating interval (Enholm et al. (2010); Parbinary strategies to scale the sonication power and/or induced side effects, making it highly appealing for clinof the heating interval et (2010); Parnot toxicity in healthy tissueappealing and therates thereby induced side effects, making it highly for clinical increase application to increase treatment success and length binary strategies to basic scale the(Enholm sonication power and/or length of the(2012)), heating interval (Enholm et al. al. (2010); Partanen et al. proportional-integral-derivative ical application to increase treatment success rates and tanen et al. (2012)), basic proportional-integral-derivative induced side effects, making itpatients highly in appealing for clinimprove quality of life for the therapy. ical application to increase treatment success rates and length ofal. the(2012)), heatingbasic interval (Enholm al. (2010); Par(PID)-based methods (Mougenot et al.et(2009)) or some tanen et proportional-integral-derivative improve quality of life for the patients in therapy. methodsbasic (Mougenot et al. (2009)) or some ical application to life increase rates and (PID)-based improve quality of for thetreatment patients insuccess therapy. tanen al. (2012)), proportional-integral-derivative hybridetform ofmethods PID and(Mougenot bang-bang (Staruch al. (PID)-based etcontrol al. (2009)) or et some A particularly interesting technology for hyperthermia is form ofmethods PID and(Mougenot bang-bangetcontrol (Staruch et al. improve qualityinteresting of life for the patientsfor in therapy. A particularly technology hyperthermia is hybrid (PID)-based (2009)) some (2011); form Lin of et PID al. (1990); Bing et control al.al.(2015)). Aor major hybrid and bang-bang (Staruch et al. magnetic-resonance-guided high-intensity focused ultraA particularly interesting technology for hyperthermia is Lin et al. (1990); Bing et al. (2015)). A major magnetic-resonance-guided high-intensity focused ultra- (2011); form of and bang-bang (Staruch ettake al. (2011); Lin et PID al.designs, (1990); Bing etiscontrol al. (2015)). A to major drawback of such however, their inability A particularly interesting technology hyperthermia is hybrid sound (MR-HIFU), which enables heatfor delivery to internal magnetic-resonance-guided high-intensity focused ultradrawback of such designs, however, isal. their inability to take sound (MR-HIFU), which enables heat delivery to internal (2011); Lin et al. (1990); Bing et (2015)). A major the body’s future thermal behavior and the restrictive acdrawback of such designs, however, is their inability to take magnetic-resonance-guided high-intensity focused ultratissues(MR-HIFU), with spatial accuracy in theheat millimeter range, com- the body’s future thermal behavior and the restrictive acsound which enables delivery to internal tissues with spatial in millimeter range, comoffuture such designs, however, isand their inability take tuator constraints into account when computing theto input, the body’s thermal behavior the restrictive acsound (MR-HIFU), which enables heat toresponse internal tissueswith with real-time spatial accuracy accuracy in the theof millimeter range, com- drawback bined monitoring thedelivery thermal tuator constraints into account whenand computing the input, bined with real-time monitoring of the thermal response body’s future affects thermal behavior theand restrictive acwhich negatively treatment duration. tuator constraints into account whenquality computing the input, tissues with real-time spatial accuracy in theofmillimeter range, com- the bined with monitoring the thermal response which treatment and ⋆ This research tuator constraintsaffects into account whenquality computing the input, which negatively negatively affects treatment quality and duration. duration. has been made possible of by the Cancer Society bined with real-time monitoring theDutch thermal response ⋆ This research has been made possible by the Dutch Cancer Society In thisnegatively respect, model predictive quality control (MPC) offers ⋆ This which affects treatment duration. and theresearch Netherlands Organisation for by Scientific Research has been made possible the Dutch Cancer (NWO) Society In this respect, model predictive control and (MPC) offers and the Netherlands Organisation for Scientific Research (NWO) superior potential for hyperthermia. This has been recIn this respect, model predictive control (MPC) offers ⋆ as This part of theirhasjoint Partnership Programme: “Technology for been made possible by the Dutch Cancer (NWO) Society and theresearch Netherlands Organisation for Scientific Research superior potential for hyperthermia. This has been recas part of their joint Partnership Programme: “Technology for In this by respect, model predictive control (MPC) superior potential for hyperthermia. Thisresults has been recognized some researchers, of which the areoffers pubOncology”. This is partially financed by the PPP Allowance as part theirproject jointOrganisation Partnership Programme: “Technology for and the of Netherlands for Scientific Research (NWO) ognized by some researchers, of which the results are pubOncology”. This project is partially financed by the PPP Allowance superior potential for hyperthermia. This has been reclished in, for example, Arora et al. (2007) for the control ognized by some researchers, of which the results are pubmade available by Top Sector Life Sciences & Health. This work is Oncology”. is partially financed by the“Technology PPP Allowance as part of This theirproject joint Partnership Programme: for lished in, for example, Arora et al. (2007) for the control of of made available by Top Sector Life Sciences & Health. This work is ognized researchers, of al. which thefor results are pubalso partially funded by Sector the European Union & via IPaCT Project. lished in,byforsome example, Arora et (2007) the control of made available byproject Top Life Sciences This work is Oncology”. This is partially financed byHealth. the PPP Allowance also partially funded by the European Union via the IPaCT Project. lished in, for example, Arora et al. (2007) for the control of also partially funded by Sector the European Union & viaHealth. the IPaCT made available by Top Life Sciences This Project. work is also partially funded by the European Union via the IPaCT Project. 2405-8963 © 2018 2018, IFAC (International Federation of Automatic Control) Copyright © IFAC 223 Hosting by Elsevier Ltd. All rights reserved. Copyright © under 2018 IFAC 223 Control. Peer review responsibility of International Federation of Automatic Copyright © 2018 IFAC 223 10.1016/j.ifacol.2018.11.012 Copyright © 2018 IFAC 223
2018 IFAC NMPC 192 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
thermal dose in a single point based on a one-dimensional temperature distribution, De Bever et al. (2014) for thermal dose control in two-dimensional systems using heavily simplified models and a fixed sonication trajectory, and more recently in Sebeke et al. (2017) for temperature control in two-dimensional systems with online adaption of the heating location and power. Although the initial results of online optimization-based control in the aforementioned works look promising, it is also clear that the full potential of MPC is far from reached. More sophisticated MPC setups specifically tailored to hyperthermia can fulfill the envisioned potential, resulting in more reliable, reproducible and highperforming 3D temperature control, which in turn may serve as an enabler for widespread and successful clinical use of MR-HIFU for hyperthermia. It is the objective of this paper to provide a significant step forward in the development of this societally highly relevant technology. To this end, an offset-free MPC design is developed for planar heating using an MR-HIFU system. A Luenbergertype observer is employed to improve the MR thermometry measurements and to construct a disturbance estimate. The latter enables the controller to compensate for the effects of constant (or slowly varying) model errors, allowing for enhanced control performance in the clinic, where such model discrepancies are inevitable. Simulations are performed to demonstrate the developed model-based algorithm’s potential in terms of achieving desired steadystate behavior despite severe plant-model mismatch. The remainder of this paper is organized as follows. In Section 2 relevant aspects related to the clinical application of mild hyperthermia are discussed in more detail to clarify the challenges involved with feedback control for this application, and to motivate the desire to develop MPC for MR-HIFU. Next, Section 3 explains how the thermal model is obtained and Section 4 covers the controller design procedure. Simulations are performed to evaluate the algorithm’s ability of removing the steady-state offset otherwise induced by model mismatch, the results of which are analyzed in Section 5. Finally, Section 6 summarizes the key achievements and corresponding observations. 2. PROBLEM DESCRIPTION This section introduces some practical aspects concerning the hyperthermia treatment and setup, and discusses the motivation for using MPC in such an application. 2.1 MR-HIFU hyperthermia setup Although the MPC setup we propose is generic in nature, in this paper it is designed for the Philips Sonalleve shown in Fig. 1, which is a commercial MR-HIFU system already being used in clinics to treat uterine fibroids and for bone pain palliation therapy. It is a combined setup of a standard MRI scanner used for noninvasive thermography and a dedicated trolley-tabletop in which an MR-compatible HIFU transducer is integrated. 2.2 HIFU applicator This system uses a phased-array HIFU transducer to generate the ultrasound waves. It consists of 256 elements 224
Fig. 1. Philips Sonalleve MR-HIFU therapy platform. focal plane focal point
transducer driving signals
tumor
patient table transducer axis
Fig. 2. Schematic of HIFU beam into focal plane in tumor area, with focal point by electronic beam steering. of which the phases and amplitudes can individually be chosen such that by interference a focal spot is created. By modulating these settings, which is referred to as electronic beam steering, the lateral location of the focal spot can be chosen anywhere on the focal plane within an 18 mm diameter circle around the transducer axis, as depicted in Fig. 2. Since the size of the focal spot is typically smaller than the tumor, it must be scanned through the treatment volume using beam steering or by moving the transducer. In this paper, we focus on beam steering in the focal plane using a predetermined set of discrete possible focus locations. 2.3 Motivation for offset-free MPC in hyperthermia To fully benefit from the hyperthermia-induced effects such as increased blood flow and oxygenation, it is crucial to maintain temperatures above 41 ◦ C in the region of interest (ROI), typically the tumor and possibly some adjacent tissue, during the entire 90 minute treatment. On the other hand, temperatures above 43 ◦ C are undesired since they may cause reverse effects such as vascular shutdown (Hijnen et al. (2014)). Furthermore, temperature elevation outside the ROI must be avoided to prevent damage to healthy tissue. MPC is considered superior to PID-based strategies for application in MR-HIFU hyperthermia, as it is able to reduce treatment time by exploiting beneficial (future) behavior, e.g., heat diffusion via conduction, and can explicitly incorporate actuator constraints such as the inability to actively remove heat from inside the body using HIFU. Unfortunately, model-based strategies are inherently accompanied by the possibility of modeling errors, which in the thermal models for hyperthermia treatment are highly likely to occur due to the high variability of spatially distributed and time/temperaturevarying tissue parameters. Inadequately accounting for such mismatches may result in insufficient heating of the target region or in the overheating of healthy tissue, thereby significantly deteriorating treatment quality.
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
acoustic properties of the intermediate tissue. In this paper, however, F is modeled by a radially symmetric twodimensional Gaussian distribution centered around the focus location rf with standard deviation σf = 2.4 mm �r − rf (t)�2 α , (3) exp − F (r, t) = 2πσf2 2σf2
MR-HIFU setup and patient
Desired acoustic power
with r ∈ Ω, t ∈ R≥0 , and scaling factor α ∈ R>0 .
MRT data MPC
193
3.2 Discrete-time state-space model
Model-based observer T setpoint, constraints
Fig. 3. Block scheme of envisioned fully automated MRHIFU feedback loop for patient treatment. Given this situation, we propose an MPC setup novel for hyperthermia inspired by an offset-free MPC algorithm as in Pannocchia and Rawlings (2003). Fig. 3 depicts the resulting feedback scheme, where the observer constructs the model-based temperature and disturbance estimates from the noninvasively obtained MR thermometry data, such that the MPC scheme can compute the optimal power distribution to be generated by the HIFU actuator. 3. THERMAL PROCESS MODELING In this section, a model of the body’s thermal response will be discussed, from which the state-space models used by the observer and MPC algorithm will be derived. 3.1 Bioheat model The tissue’s temperature change in response to the applied heating is modeled using Pennes’ bioheat equation, see Pennes (1948), given by ∂T (r, t) = k∇2 T (r, t) + Q(r, t) − wb (r)cb (T (r, t) − Tb ), ρc ∂t (1) where T : Ω × R≥0 → R is the temperature profile, such that T (r, t) denotes the temperature at time t ∈ R≥0 and location r = [rx , ry ]⊤ ∈ Ω ⊂ R2 with Ω being the patient domain in the focal plane. The parameters ρ, c, k ∈ R>0 denote the tissue’s volumetric mass density, specific heat capacity, and thermal conductivity, and wb : Ω → R≥0 and cb ∈ R>0 the blood’s perfusion coefficient and specific heat capacity, respectively. In this paper, the tissue is assumed to be homogeneous with respect to most parameters (ρ, k, c, and cb ), but not the blood perfusion rate wb (r), which is assumed significant only in tissue in the direct vicinity of blood vessels, and zero elsewhere. Note, however, that our MPC setup can be extended to incorporate the fully inhomogeneous case as well. The power deposition density described by Q : Ω × R≥0 → R≥0 depends on the acoustic deposition intensity profile F : Ω × R≥0 → R≥0 and linearly scales with the sonication power P : R≥0 → R≥0 according to Q(r, t) = F (r, t)P (t), r ∈ Ω, t ∈ R≥0 . (2) In reality, for a given focus location rf (t) at some time t, the intensity F (r, t) also depends on factors such as the 225
Assuming zero out-of-plane interaction, (1)-(3) are discretized spatially using the central difference scheme on a 39 × 39 grid with a voxel size of 2 × 2 mm2 and temporally using forward Euler with a sampling interval of Ts = 2.5 s, resulting in the discrete-time state-space model xk+1 = Ap xk + Bp uk , (4a) (4b) yk = xk + vk , with k ∈ N connecting to real time tk = kTs and the states xk ∈ Rn , n = 392 = 1521 representing the voxel temperatures. We select 48 voxels within the ROI, at the center of which we allow sonication by steering the focal point rf to this location. These locations are referred to as sonication points. The input uk ∈ Rm with m = 48 is the applied sonication power averaged over the sampling interval at each of these sonication points. Each individual column of Bp ∈ Rn×m captures the spatially discretized power deposition as in (2) and (3) corresponding to rf at a single sonication point, and its temporally discretized relation to voxel temperature change according to (1). The spatial discretization is chosen in such a manner that the centers of the voxels coincide with the points measured by MR thermometry. Also, the sampling interval Ts is larger than the time required for measurement, implying that at each discrete time k ∈ N a full new temperature map is available and thus the output yk ∈ Rn can indeed be modeled as a full (noisy) state measurement as in (4b). In practice, the MR data is corrupted by noise, which for the purpose of this paper can be approximated by uncorrelated zero-mean Gaussian white vk ∈ Rn , √ noise ◦ k ∈ N, with standard deviation σv = 0.2 C. 4. CONTROLLER DESIGN In this section, the state-space models used by the controller and observer are given first. Then, the control objectives as described in Subsection 2.3 are formulated as a constrained optimization problem. 4.1 Prediction model To enable the controller to compensate for the steadystate offset induced by model mismatch or slowly varying disturbances, integral action is incorporated using disturbance estimation as described in Pannocchia and Rawlings (2003). The state-space representation of the prediction model is given by xi+1|k = Axi|k + Bui|k + di|k , (5a) di+1|k = di|k , (5b) yi|k = xi|k , (5c) where xi|k , ui|k , yi|k ∈ Rn denote the predicted states, inputs, and outputs, respectively, at i ∈ N time steps ahead
2018 IFAC NMPC 194 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
of the prediction sequence’s starting time k. Furthermore, di|k ∈ Rn denotes the virtual state disturbance, which is assumed to be constant over the prediction horizon. Note that we allow model mismatch, i.e., A �= Ap and B �= Bp , and that we did not incorporate measurement noise in the prediction model since E(vk ) = 0. As proposed in Theorem 1 of Pannocchia and Rawlings (2003), we first verify that the nominal model (i.e., the model without additional disturbance di|k ) is detectable (in fact, observable). Second, we evaluate the rank condition stated in the theorem, which using our disturbance model becomes I − A −I rank = 2n, (6) I 0
to conclude that (5) is detectable, and therefore an asymptotically stable estimator exists given the proposed disturbance model. This property also guarantees that for a stable closed-loop system, a feasible steady-state setpoint, and constant disturbances, there will be zero offset between the realized and desired temperature distributions in case no constraints are active at steady state. 4.2 State and disturbance estimation To facilitate this offset-free implementation, and to simultaneously obtain improved temperature estimates with respect to the noise-corrupted MR thermometry readings, we design the state and disturbance estimator x ˆk = Aˆ xk−1 + Buk−1 + dˆk−1 + Lx (yk − yˆ− ), (7a) k
(7b) dˆk = dˆk−1 + Ld (yk − yˆk− ), where yˆk− = Aˆ xk−1 + Buk−1 + dˆk−1 (8) denotes the model-based measurement estimate at time k before applying the correction by output injection to obtain x ˆk . Note that we choose the observer model (7) equal to the prediction model (5). Since we have full (noisy) state measurement, the system (5) is observable, allowing for arbitrary pole placement in (7) via Lx and Ld . We choose the diagonal matrices 1 1 In , Lx = In , Ld = (9) 3 30 which cause the state and disturbance estimation errors to converge to within 5% of their initial value in approximately 22.5 and 225 seconds, respectively. 4.3 Optimization problem
T [◦ C] 43 42 41 40
T +ǫ T −ǫ
Transducer axis ǫ
ǫ
37
R
r2Ω
18 mm
S 32 mm
Fig. 4. Schematic cross-section of the temperature objectives for the circular regions R and S centered at the transducer axis. The maximum violations ǫ and ǫ are shown for some over- and underheated temperature distributions such that T1 ≤ T + ǫ and T2 ≥ T − ǫ. lower value on Ω \ S to safeguard healthy tissue. Within the aforementioned temperature range, preference is given to a flat temperature distribution in the ROI, which is described by the reference temperature Tr . To formulate these objectives mathematically, let us denote the elements of x corresponding to a point within R by xR = Hx ∈ RnR , where H ∈ {0, 1}nR×n is a matrix with exactly one 1 in each row. Furthermore, let xr ∈ RnR and x, x ∈ Rn denote the voxel-wise temperature reference and lower/upper bounds corresponding to the values of Tr , T , and T on the discrete voxel locations, respectively. The voxels’ temperature deviations with respect to the reference in the ROI are then given by xR − xr . The maximum violations of the temperature bounds follow from ǫ = �max(x − x, 0)�∞ ≥ 0 and ǫ = �max(x − x, 0)�∞ ≥ 0, where the maximum operator is used component-wise, and are collected in ǫ = [ǫ, ǫ]⊤ ∈ R2≥0 . The control objective is now formulated as the constrained optimization problem N −1 min ℓ(xi|k , ui|k , ǫi|k ) + ℓ(xN |k , 0, ǫN |k ) (10a) u , i|k
i∈[0,N −1]
i=0
with, for some i, k ∈ N, the stage cost r r ⊤ R ℓ(xi|k , ui|k , ǫi|k ) = (xR i|k − x ) Q(xi|k − x ) ⊤ + u⊤ i|k Rui|k + fǫ ǫi|k ,
The temperature objectives have been discussed in Subsection 2.3 and are schematically depicted in Fig. 4 in cross-section perspective. T : Ω → R and T : Ω → R represent the (location-dependent) lower and upper temperature bounds, respectively, defining the desired temperature range (green) over the entire patient domain Ω. Underand overheating occurs for temperatures within the (light) blue and (light) red areas in the figure, respectively. The concentric circular regions R and S denote the ROI and the region additionally including the tissue surrounding the ROI, i.e., R ⊂ S ⊂ Ω. T is nonzero only within R where sufficient heating is desired . T features an elevated plateau on S, preventing reversal of the temperaturedependent mechanisms due to overheating, and has a 226
Tr T;T T1 T2
subject to x − 1ǫi|k xi+1|k di+1|k ui|k
(10b)
i ∈ [0, N ], (10c) ≤ xi|k ≤ x + 1ǫi|k , = Axi|k + Bui|k + di|k , i ∈ [0, N ], (10d) = di|k , i ∈ [0, N − 1], (10e) ≥ 0, i ∈ [0, N − 1], (10f)
1⊤ ui|k ≤ utot ,
i ∈ [0, N − 1]. (10g)
Inequality (10c), in which 1 denotes an all-ones vector of suitable dimension, captures the temperature bounds as soft constraints (note that ǫi|k = [ǫi|k , ǫi|k ]⊤ ), (10d)(10e) ensure satisfaction of the dynamics according to (5) with the prediction sequence’s initial conditions given by ˆk and d0|k = dˆk , and finally (10f)-(10g) describe x0|k = x
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
which are normalized with respect to the number of corresponding variables for more intuitive balancing of the objectives’ relative contribution to the cost function. The critical importance of achieving temperatures within the desired range is expressed by the linear term (which is possible since ǫ ≥ 0 by definition) with high weighting. Using only the maximum violation, as opposed to penalizing each voxel individually, improves temperature homogeneity during the initial heat-up phase. To also encourage a flat temperature distribution (within the temperature bounds) in the ROI at steady steady, a quadratic term with relatively lower cost is added to penalize this tracking error. In addition to the aforementioned control objectives, a small weight on the input power is incorporated to favor solutions that involve less heating.
195
43 Temp. in ROI [◦ C]
the actuator constraints with utot = 60 W being the maximum allowable power being applied to the treated region. Furthermore, a horizon of N = 5 is used, and the weights are chosen to be 0.01 10 100 I, fǫ = I, R = , (11) Q= 200 nR m
42 41 xR avg R xR min , xmax R x ˆavg x ˆR ˆR max min , x
40 39 38 37
0
100
200
300 Time [s]
400
500
600
400
500
600
(a)
43 Temp. in ROI [◦ C]
2018 IFAC NMPC Madison, WI, USA, August 19-22, 2018
42 41 40 39 38 37
0
100
200
5. PERFORMANCE ANALYSIS
300 Time [s] (b)
Next, we perform a simulation in which using the nominal model results in temperature overestimation. To this end, the thermal conductivity k in (1) is intentionally underestimated and the power deposition intensity’s scaling factor α in (3) is overestimated in the prediction and estimation models, both by a factor of 5, affecting A and B, respectively, in (5) and (7). Additionally, the observer and controller models are set to assume a homogeneous perfusion rate wb (r) = wb , thereby neglecting the local influence of blood vessels present in the plant. The results are shown in Fig. 5b. Initially, the disturbance estimator is disabled, i.e., dˆk = 0 and Ld = 0 for k ∈ N corresponding to tk < 300 s. This causes a mismatch between x ˆR and xR , resulting in a downward offset in the realized xR with respect to its reference xr . From tk = 300 s (vertical dashed black line) onwards, the offset-free algorithm is active by setting Ld as ˆ the achieved in (9). At tk = 600 s, after convergence of d, temperature distribution coincides with that of Fig. 5a, illustrating how offset-free MPC can significantly improve hyperthermia treatment quality under model error. This is more clearly illustrated in Fig. 6, where similarly to Fig. 4 cross-sections of the temperature objectives and distributions through the center of R (at rx = 0.128 m) 227
Fig. 5. The average (solid) and lowest/highest (dashed) voxel temperatures in the ROI of the plant xR (black) and observer x ˆR (gray) in case of (a) no plant-model mismatch, and (b) significant model error where disturbance estimation is enabled at t = 300 s. Temperature [◦ C]
A simulation study was carried out to verify the effectiveness of the proposed control algorithm in terms of compensating for the effects of plant-model mismatch for the MR-HIFU system under consideration. For reference, a simulation without model error is performed first, i.e., where A = Ap and B = Bp , of which the results are shown in Fig. 5a. Initially, the entire patient volume is assumed to be 37 ◦ C. The estimated temperature x ˆR (gray) closely R resembles its true counterpart x (black) in the ROI, as is shown in the figure by their lowest/highest (dashed) and average value (solid), enabling the controller to accurately heat the target region in a satisfactory manner. The current MPC’s worst-case computation time is found to be 1 s, indicating it is sufficiently fast given Ts = 2.5 s.
43 42 41 40 39 38 37 0.09 0.10
(a) (b) (c) Fi
0.11
0.12
0.13 ry [m]
0.14
0.15
0.16
Fig. 6. Voxel temperature xk at cross-section through center of ROI without plant-model mismatch at tk = 600 s (a), and with mismatch at tk = 300 s (b) and tk = 600 s (c), alongside the acoustic intensity around in-plane sonication points Fi . The gray areas indicate the cross-section voxels with increased perfusion. are shown, alongside the (scaled) acoustic intensities Fi = F (rf,i ), i = 1, . . . , 6, from (3) with center locations rf,i at the six in-plane sonication points. The temperature profile corresponding to Fig. 5a at tk = 600 s (black) shows the MPC’s optimal performance without model error. It is within the desired temperature range, i.e., x ≤ x ≤ x, and tends to xr in R, however due to the system restrictions it is not a perfectly flat distribution. That is, the acoustic intensity profiles (blue) do not admit further heating of the areas with locally increased perfusion (gray) without causing additional overheating outside these areas. As previously found in Fig. 5b at tk = 300 s, Fig. 6 shows that while introducing plantmodel mismatch with disturbance estimation disabled
ry [m]
2018 IFAC NMPC 196 Madison, WI, USA, August 19-22, 2018
D.A. Deenen et al. / IFAC PapersOnLine 51-20 (2018) 191–196
0.16
0.6
0.15
0.4
0.14
0.2
0.13
0
0.12
−0.2
0.11 0.10
−0.4
0.09
−0.6 0.10
0.14 0.12 rx [m]
0.16
Fig. 7. Disturbance estimate dˆ at time tk = 600 s, with additionally the edges of the circular R (solid black) and S (dashed black), sonication points (black dots), and voxels with increased perfusion (red squares). causes insufficient heating (solid red), optimal performance is recovered by the offset-free algorithm after convergence at tk = 600 s (dashed red). The disturbance estimate after convergence at tk = 600 s is visualized in Fig. 7. Indeed, overestimation occurs in the entire ROI due to the mismatch in k and α, which ˆ Moreover, the is compensated for by negative values in d. figure shows that the disturbance estimator correctly identifies the voxels in which the locally increased perfusion is neglected (red outline), and thus is able to counteract the corresponding effects. In the figure, some voxels outside S have a slightly positive disturbance estimate, which is a result of the propagation of measurement noise into dˆ via the output injection in (7). 6. CONCLUSION In this paper, an offset-free MPC scheme is presented to enhance control performance of an MR-HIFU hyperthermia system used in cancer treatment. The proposed design achieves more accurate 2D temperature control than the techniques currently used in clinics. Specifically, the model-based algorithm is able to explicitly take into account (a) the restrictive input constraints and (b) the body’s future thermal response, thereby enabling more optimal heating and shorter treatment times. Using disturbance estimation, the control design eliminates the steady-state offset otherwise resulting from modeling error, thus enhancing hyperthermia treatment quality. The MPC setup’s effectiveness in terms of compensating for unmodeled perfusion and inaccurate conduction/heating parameters is demonstrated using simulation studies. REFERENCES Arora, D., Skliar, M., Cooley, D., and Roemer, R.B. (2007). Constrained Predictive Control of Thermal Therapies for Minimum-Time Delivery of Thermal Dose. IEEE Trans. Control Syst. Technol., 15(6), 1030– 1037. Bing, C., Nofiele, J., Staruch, R., et al. (2015). Localised hyperthermia in rodent models using an MRIcompatible high-intensity focused ultrasound system. Int. J. Hyperth., 31(8), 813–822. 228
Datta, N., Ord´on ˜ ez, S.G., Gaipl, U., et al. (2015). Local hyperthermia combined with radiotherapy and-/or chemotherapy: Recent advances and promises for the future. Cancer Treat. Rev., 41(9), 742–753. De Bever, J., Todd, N., Payne, A., Christensen, D.A., and Roemer, R.B. (2014). Adaptive model-predictive controller for magnetic resonance guided focused ultrasound therapy. Int. J. Hyperth., 30(7), 456–470. Enholm, J., Kohler, M., Quesson, B., et al. (2010). Improved Volumetric MR-HIFU Ablation by Robust Binary Feedback Control. IEEE Trans. Biomed. Eng., 57(1), 103–113. Franckena, M., Fatehi, D., De Bruijne, M., et al. (2009). Hyperthermia dose-effect relationship in 420 patients with cervical cancer treated with combined radiotherapy and hyperthermia. Eur. J. Cancer, 45(11), 1969–1978. Hijnen, N., Langereis, S., and Gr¨ ull, H. (2014). Magnetic resonance guided high-intensity focused ultrasound for image-guided temperature-induced drug delivery. Adv. Drug Deliv. Rev., 72, 65–81. Issels, R.D., Lindner, L.H., Verweij, J., et al. (2018). Effect of Neoadjuvant Chemotherapy Plus Regional Hyperthermia on Long-term Outcomes Among Patients With Localized High-Risk Soft Tissue Sarcoma. JAMA Oncol., 4(4), 483–492. Lin, W.L., Roemer, R.B., and Hynynen, K. (1990). Theoretical and experimental evaluation of a temperature controller for scanned focused ultrasound hyperthermia. Med. Phys., 17(4), 615–625. Mallory, M., Gogineni, E., Jones, G.C., Greer, L., and Simone, C.B. (2016). Therapeutic hyperthermia: The old, the new, and the upcoming. Crit. Rev. Oncol. Hematol., 97, 56–64. Maloney, E. and Hwang, J.H. (2015). Emerging HIFU applications in cancer therapy. Int. J. Hyperth., 31(3), 302–309. Mougenot, C., Quesson, B., De Senneville, B.D., et al. (2009). Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU). Magn. Reson. Med., 61(3), 603–614. Pannocchia, G. and Rawlings, J.B. (2003). Disturbance models for offset-free model-predictive control. AIChE J., 49(2), 426–437. Partanen, A., Yarmolenko, P.S., Viitala, A., et al. (2012). Mild hyperthermia with magnetic resonance-guided high-intensity focused ultrasound for applications in drug delivery. Int. J. Hyperth., 28(4), 320–336. Pennes, H. (1948). Analysis of Tissue and Arterial Blood Temperatures in the Resting Human Forearm. J. Appl. Physiol., 1(2), 93–122. Sebeke, L., Luo, X., Heijman, E., et al. (2017). Predictionbased controller for MR-HIFU mediated hyperthermia. In Eur. Symp. Focus. Ultrasound Ther., 23. Staruch, R., Chopra, R., and Hynynen, K. (2011). Localised drug release using MRI-controlled focused ultrasound hyperthermia. Int. J. Hyperth., 27(2), 156–171. Staruch, R.M., Hynynen, K., and Chopra, R. (2015). Hyperthermia-mediated doxorubicin release from thermosensitive liposomes using MR-HIFU: Therapeutic effect in rabbit Vx2 tumours. Int. J. Hyperth., 31(2), 118–133.