Development of robust/predictive control strategies for image-guided ablative treatments using a minimally invasive ultrasound applicator

Abstract

Purpose: One important challenge in image-guided ablative therapies is the effect of heat diffusion which can cause damage to surrounding organs and limit the ability to achieve a conformal pattern of thermal damage. Furthermore, tissue properties such as perfusion and energy absorption can be dynamic and difficult to measure. This paper attempts to address these problems by proposing new control methods. Materials and methods: A novel predictive approach was developed to compensate for the effect of heat diffusion using a minimally invasive rotating ultrasound heating applicator for ablative therapy. This method can be merged into any closed-loop control strategy. A binary controller, a previously developed adaptive proportional-integral controller, and a model reference adaptive controller were employed and compared, all with the predictive element incorporated. The reason for choosing these controllers was that none of them needed a model of the tissue or exact values of their parameters. Results: The effectiveness of these controllers was demonstrated through both simulation and experimental studies. The results were consistent and demonstrated equivalent performance between controllers. Conclusions: The dominant influence on radial targeting accuracy was the prediction element described in this paper. A binary controller with a predictive element may provide the best balance of performance and simplicity for this application.

• Heat diffusion
• magnetic resonance imaging (MRI)
• minimally invasive
• prediction
• robust temperature control
• ultrasound interstitial thermal therapy

Introduction

Image-guided ultrasound therapy is a powerful method of cancer therapy in which high intensity ultrasound energy is used to coagulate (or cook) a target region of tumour from outside the body or within body cavities [1,2]. This technology can be used in soft tissue organ sites where a target volume of tissue is amenable to thermal coagulation. Across all of these organs there are common aspects that the technology must achieve with respect to accuracy of spatial heating, but also unique characteristics with respect to safety of surrounding structures. This requires accurate control over spatial and temporal deposition of energy to regulate the temperature. To deliver ultrasound energy to the targeted tumour there are two different technologies: (1) external high-intensity focused ultrasound (HIFU) therapy and (2) interstitial ultrasound therapy. In external HIFU therapy the ultrasound transducer is placed outside the body and energy is transmitted into the tumour. In minimally invasive ultrasound therapy the transducer is inserted directly into the tumour or an adjacent body cavity. External HIFU is often difficult for soft tissue targets located behind bony structures or gas, and in these situations, minimally invasive approaches can be more suitable. Furthermore, these approaches are often more efficient for volumetric tissue ablation due to the direct deposition of ultrasound energy in tissue and the lack of need to protect intervening tissues.

Feedback control techniques can be used to regulate the temperature and have many advantages over open-loop approaches because they can compensate for unpredictable disturbances and time-varying changes in tissue properties such as perfusion and ultrasound absorption. To monitor the temperature in real time for the feedback control loop, magnetic resonance (MR) thermometry using the proton resonant frequency shift (PRFS) method can be used [3].

One of the problems in MR-guided thermal therapy is the effect of heat diffusion which can cause damage to surrounding structures. Thermal dose in tissue can increase even when ultrasound power is turned off. A model predictive control technique was used to overcome this problem [4]. The method is limited to controlling the temperature at a single point using an externally focused ultrasound source. Furthermore, the method is highly dependent on the model of the tissue. The primary objective of this paper was to propose a novel approach for solving the problem of heat diffusion for interstitial ultrasound therapy while trying to control the temperature in a 3D treatment volume. The proposed algorithm is very general and can be incorporated into any feedback control strategy, and is not dependent on the model of the tissue.

Another main difficulty is the uncertainty in the parameters used for tuning the fixed gain controllers such as the ultrasound absorption or blood perfusion [5]. These properties are known to be spatially heterogeneous and difficult to measure in tissue. Furthermore, there is no theoretically rigorous and experimentally validated bioheat transfer model [6]. The second objective of this paper was to find an effective approach to control the temperature in a dynamic environment. Having a controller which does not depend on the model of the system can be beneficial for this application.

In this context, temperature control using adaptive control techniques such as a self-tuning regulator (STR) and model reference adaptive control (MRAC) have been reported by Sun et al. [7]. In Hey et al. [8], an adaptive proportional integral derivative (PID) controller was developed to control the temperature. However, these applications were only limited to a single point temperature control and used an externally focused ultrasound transducer. Also, the effect of heat diffusion was not considered in these studies. The third objective of this paper was to investigate the effectiveness of different controller strategies towards 3D thermal therapy considering the effect of heat diffusion on controllers’ performance.

Although the results presented in this paper can be used for other interstitial treatments using a rotating ultrasound beam in brain or liver tissue, the focus in this study is on a transurethral treatment for prostate cancer. The experimental set-up developed within our group for MRI-controlled transurethral ultrasound therapy of prostate cancer consists of a multi-element ultrasound heating applicator operated under rotational control within standard MRI [9,10,11]. In similar work by other groups [12,13], the power was set to a constant value and the rotation rate was changed manually. In our group, the rotation rate of the applicator and the frequency and power output from each element on the transducer were automatically adjusted such that a desired temperature value was achieved along the pre-defined target boundary. This approach results in a multi-input single-output (MISO) control problem, making it more challenging [10].

Three different controllers are selected for a comparison study. A MRAC, an adaptive PI controller [11] and a binary controller are developed and compared; all with the predictive element incorporated. The reason for choosing these controllers was that none of them needed a model of the tissue or exact values of its parameters. In this novel predictive approach we not only considered the boundary temperature along the beam path, but also incorporated into the error signal previous weighted temperature values over a fixed window. The effectiveness of the aforementioned controllers employing the prediction element are validated and compared through both simulation and experimental studies.

Experimental set-up

Heating experiments were performed using a MRI-compatible treatment delivery system developed for research within our group, comprising a multi-element ultrasound heating applicator, a rotational piezoelectric motor and radio-frequency (RF) electronics, as described in detail by Chopra et al. [9,10] (see Figure 1).

Figure 1. (A) The concept of MRI-controlled transurethral ultrasound therapy with multiple collimated high-intensity ultrasound beams. (B) Photograph of an MRI-compatible transurethral ultrasound applicator used to deliver high-intensity ultrasound energy to the prostate gland. (C) Photograph of the experimental system used in this study attached to the patient table of a clinical 3-T MR imager.

MR thermometry was performed using the PRFS method [3] continuously during ultrasound delivery and rotation to measure the spatial temperature distribution in the prostate gland and surrounding tissues. This concept is shown in Figure 1(A). Multiple MR images were acquired in the planes of the transducers (dotted lines) continuously during treatment. The temperature maps were analysed by the treatment delivery system as soon as they were acquired (every 5 s). This feature of simultaneous motion and imaging was a central aspect of the feedback control for this technology. Based on the temperatures measured along the direction of the ultrasound beam, the rotation rate of the applicator, the frequency, and the power output for each element on the transducer were adjusted so that a desired temperature threshold was achieved along the pre-defined target boundary.

Ultrasound and thermal calculations

This study used computer simulations, implemented in C++, to evaluate the accuracy with which volumetric regions of thermal coagulation could be shaped to target boundaries using multi-element transurethral ultrasound heating applicators and feedback control. In order to evaluate the performance of the control algorithms on realistic target boundaries, 3D models were created by segmenting MR images of prostate cancer patients. The computational domain was a large 16 × 16 × 16-cm³ volume.

The complex 3D acoustic pressure distribution from a single planar rectangular transducer element was calculated for each frequency using an approximate solution to the Rayleigh–Sommerfeld integral by Ocheltree and Frizzell [14]. Tissue temperature dynamics due to ultrasound power deposition, heat diffusion, and thermal homeostasis from blood perfusion were modelled using an explicit 3D finite difference time-domain solution to Pennes’ bioheat transfer equation [15]. The equation is given by where k_t is the thermal conductivity of tissue (W m⁻¹° C⁻¹), is the temperature distribution in tissue, T_b is the temperature of blood, Q(r,t) is the deposited acoustic power (W m⁻³), w_b is the blood perfusion (kg m⁻³ s⁻¹), c_b is the specific heat capacity of blood (J kg⁻¹° C⁻¹), is the specific heat capacity of tissue and ρ is the density of tissue (kg m⁻³).

The values for tissue parameters used in the simulations are shown in Table 1 [10].

Table 1. Physical parameters used in acoustic calculation and bio-thermal simulation [10].

Parameter	Value
Density, ρ	1000 kg m^−3
Specific heat capacity, cb = ct	3700 J kg^−1° C^−1
Blood perfusion, wb	5.3 kg m^−3 s^−1
Thermal conductivity, kt	0.5 W m^−1° C^−1
Attenuation, μatt	5.76 Npm^−1MHz^−1
Ultrasound absorption, αabs	7.4 Npm^−1MHz^−1
Speed of sound, c	1500 m s^−1

Description of a novel predictive approach for transurethral cancer therapy

In the transurethral ultrasound cancer therapy being developed by our group for prostate cancer, a multi-element transducer was positioned inside a continuous 3D target boundary of the prostate. Figure 2 shows this concept as well as the typical temperature profile along the direction of heating of an ultrasound transducer element.

Figure 2. Typical temperature profile along the direction of heating of an ultrasound transducer element.

The temperature rises quickly from T_u (the temperature of water flowing through the transurethral device), reaching a maximum value T_max approximately 5–10 mm from the transducer surface, then decreasing and reaching ambient body temperature within 20–30 mm (depending on power and frequency). An MR temperature map centred on each transducer element measures temperature changes generated within the target boundary in the plane of rotation of each element. The distance from each transducer element, j, to the 3D target boundary along the direction of heating was defined as the target radius, r_j, with temperature Tr_j. The value for critical temperature (T_c) was chosen to be 52 °C in this study based on preliminary studies in humans [16].

Thermal treatments can be controlled based on temperature or thermal dose. The main goal of an ablative treatment is to achieve a thermal threshold that achieves thermal coagulation. A T_c of approximately 52 °C [16] or a thermal dose of approximately 240 CEM43 are both appropriate thresholds that can be used. Therefore control of treatment can be performed with either threshold, so a controller could be developed based on either thermal dose or temperature; we chose to focus on temperature in this study. One advantage of using temperature for control is the fact that the uncertainty in measurements obtained with MR thermometry remained normally distributed. The thermal dose equation is an exponential integral relationship; therefore, the normally distributed uncertainty associated with temperature measurements becomes non-linear, resulting in an overestimation of thermal dose. This must be corrected, otherwise the erroneous thermal dose estimates can result in undertreatment.

The aim of a well-designed control algorithm is to reduce the temperature difference between T_c, the desired temperature at the boundary, and T_rj (e = T_c−T_rj) by adjusting the power (P), frequency (f) and rotation rate (w) of each transducer element. For simplicity we did not index the e and P variables throughout this paper.

If T_rj exceeds T_c, this results in overtreatment. We have observed overtreatment in simulations and experiments, even in situations where T_rj was accurately controlled to reach T_c along the beam direction. We attributed this overtreatment to the effects of heat diffusion causing slight increase in temperatures above T_c behind the beam. To overcome this source of overtreatment, a novel prediction approach was developed in this section. Consider the geometry of a typical prostate boundary depicted in Figure 3. The solid arrow is the current position of the applicator at time t and the dash arrows are assumed to be the positions of the applicator in previous times at different angles. We can collect values of the temperature in a fixed window (angle size, n) at current time t, T_rj (θ_i, t), i = 1,… ,n.

Figure 3. Geometrical description of the proposed algorithm.

In this new approach we not only considered the boundary temperature along the beam path, T_rj (θ₁, t), , but also incorporated into the error signal previous weighted temperature values over a fixed window. More specifically, instead of having the error as e = T_c −T_rj (θ₁, t), the error is augmented as where α is the weighting/forgetting factor with 0 ≤ α ≤ 1, T_c is the desired temperature value we want to reach at the boundary, and j is the number of ultrasound elements. θ_i can vary from 0° to any desired value. Note that the angular sector we look back to collect different temperatures values is fixed throughout the treatment in this initial version of the algorithm. One could consider making this window dynamic based on the anatomy of the gland, or other criteria.

Theoretical analysis

We assumed that the error under control goes to zero. However, this was not very beneficial in this application since increases due to heat diffusion. We further assumed that the amount of this increase is equal to: , i.e. .

The error equation then became where d_max is the maximum of all. We further assumed that reaches . The error equation then became

This suggests that α should be a small value. For example if α = 1 then the error is approximately equal to 0 and after heat diffusion, it will increase and cause some overtreatment. On the other hand a small value (α = 0), leads to an error close to −d_max and after heat diffusion the error roughly increases towards a zero value. In this case, the closer d_max gets to our prediction of d₁, the better treatment we would have. In order to have a better prediction for d₁ the previous data should be collected from the points which are close to the θ₁. In the next section, we further justify our claim by simulation results.

Control algorithms and their extension to 3D ultrasound cancer therapy

The aim of this paper was to find a suitable controller for interstitial ultrasound cancer therapy where an ultrasound applicator is inserted into the body. Three separate variables (power, frequency as well as the rotation rate of the ultrasound applicator) need to be controlled to regulate the temperature in this application. To control the power, three different control strategies were investigated and are briefly explained here. For simplicity a first order system model is considered as where u and y are the scalar control and output signals respectively and a,b are the unknown parameters of the system. Goharrizi et al. [10], showed that a first order system can sufficiently model the temperature rise due to ultrasound power. The error can be defined as e = y_d − y where y_d is the desired value for y. The controllers are explained next.

Binary controller

The control signal, u, for the binary controller is calculated as where u_max is the maximum value for the control signal which can be determined by the maximum power that each ultrasound element can produce without damaging it.

Adaptive proportional integral controller

For the adaptive PI controller the control signal is determined by Goharrizi et al. [11]: where K_p(t) is the proportional gain and K_I(t) is the integral gain. By selecting these adaptation laws for K_p and K_I, the stability of the closed-loop system can be guaranteed: where γ₁ > 0, γ₂ > 0 are the adaptation gains for the controller.

Model reference adaptive controller

In MRAC the desired response, y_d, is determined by a reference model. For Equation 5 the reference model of [17] is assumed: Where u_d is the desired value that the system output should reach. The control signal is as in Astrom and Wittenmark [17]:

By the Lyapunov stability theorem it can be shown that if the adaptation laws for and are chosen as follows the stability of the closed-loop system can be guaranteed: where γ > 0 is the adaptation gain.

As mentioned previously, we have a multi-input single-output control problem in this application. The output to be controlled is the temperature at the targeted boundary, T_rj, and the inputs are power (P), frequency (f) and rotation rate (ω).

First, to adjust the power, the controllers explained here were employed. The power was considered to be the input to the system (u in Equation 5) and the targeted boundary, T_rj, was considered to be the output, y, in Equation 5. The power for each element changed between zero and the maximum acoustic power, P_max, which was 4 W. The power was also turned off when the maximum temperature (T_max) exceeded a threshold value, T_h (typically 90 °C), to avoid boiling in tissue. It was also turned off if the boundary temperature exceeded the critical temperature to avoid overtreatment (e < 0). Therefore the equation of power for binary, adaptive PI and MRAC controllers turned out to be as follows, respectively:

For adaptive PI controllers and were calculated as:

The values for γ₁, γ₂ were chosen to be 0.002. For the MRAC the β₁ and β₂ were calculated as:

The value for γ was chosen to be 0.001.

Using a maximum acoustic power of 4 W, it took approximately 60 s to elevate the temperature from 37 °C to 52 °C at the beginning of treatment. In order to make sure that the reference model reached 52 °C in less time without having any overtreatment, the reference model for MRAC was chosen to be:

Another challenge in this work was to control the rotation rate. We proposed a simple and effective binary controller to adjust the rotation rate of the device and relate it to the error. Using the logic that for a high error value we needed to slow down to give the controller more time to raise the temperature, and for a small error value the applicator should speed up to avoid overheating the tissue,the following controller was proposed.

Note that the values above were found to be appropriate for fast treatment as well as accuracy. The rotation rate is limited between the maximum rotation rate, ω_max, and the minimum rotation rate ω_min. These values were set to be 40°/min and 8°/min respectively based on the finite imaging time (for the maximum rotation rate) and limitations of the rotational motor (for the lower value). Note that the above values were found to be appropriate to have a fast treatment as well as good accuracy. The frequency of the applicator was also adjusted based on the radius at which the temperature should be controlled, and is defined as

Here r_f is the target radius where the frequency was switched between the first and third harmonic in order to optimise treatment time [16]. The resonant frequencies, f_high and f_low, of the transducer were 4.5 and 14.5 MHz, based on the frequencies of existing devices, and the r_f is set to be 18 mm. r_j values are the prostate boundaries which are known beforehand.

Sensitivity analysis

In this section we address how the value of and the length of window (angle size), , can affect the results of treatment through simulation results. This approach gets feedback from the previous output of the controller treatment. If we got any overtreatment for regions close to the current position of the applicator the power was decreased, and if we got undertreatment for those regions then the power was increased. Referring to Equation 2, one can see with a small value of α more weight is given to the previous temperature values in the error signal. If the previous temperature values are more than the desired temperature value T_c the error signal decreases, which in turn decrease the power and vice versa. This also takes into account the effect of heat diffusion, cross-talk between adjacent elements, as well as changes in perfusion. The rate of heat diffusion in different locations can vary as a result of changes in tissue parameters. Since we take into account the previous output of the treatment at a location close to the current location of the applicator we can compensate for these disturbances inherited in this application. For different regions heat diffusion and radius can be changed. Due to the above reasons, for collecting the previous values of the temperature it is best not to go far away from the current position of the applicator. Furthermore, the angle size at which the previous values of the temperature should be collected is better to cover the ultrasound beam width. Throughout our experiments, we have found an angle size, , of 10–15° would be appropriate for this application. Therefore, first we gather all the temperature values at 15 degree behind the current location of rotating applicator. The mean values as well as the standard deviation of the error are calculated. The error is the difference between the prostate boundary and the boundary at which the temperature reaches the desired temperature (radial error). In this study four different patients are considered. All prostate geometries were segmented from a recent human feasibility study of patients undergoing transurethral ultrasound therapy. The data is the average over all patients. The results presented in this section are obtained using the binary controller; however, the same observation was made with other controllers. The results are shown in Figure 4. As α increased, the mean value of the error increased. This was in line with our expectations. The standard deviation of the error also decreased at first and started increasing again as α increased. To have a small mean value of error as well as less variation in error an α value of 0.3 was found to be appropriate.

Figure 4. The mean value as well as the standard deviation of the error for different values of α across all patients.

Next, to validate our hypothesis about the selection of the size of window the value of α was set to 0.3 and the influence of different sizes for the window, n, was investigated. The size of window, n, was varied from 5–40°. The mean values as well as the standard deviation of the error for different values of window size, n, were plotted in Figure 5. The mean value of the error decreased but the variation of the error increased as the window length increased. As shown, a value between 10° and 15° worked best. This validates our hypothesis that a window approximately equivalent to the width of the ultrasound beam is valid.

Figure 5. The mean value as well as the standard deviation of the error for different values of window size, n, across all patients.

The results showed the same trends for each individual patient. For example, we plotted the conformal heating pattern for one patient in Figure 6 having different values for α and the angle window of 15°. As shown, for α = 1 there was a significant overtreatment and α = 0.3 resulted in the best treatment compared to other values of α. This is in line with the results found before in this paper. The white contours in Figure 6 depict the temperature of 42 °C.

Figure 6. The pattern of conformal thermal heating for one element of ultrasound transducer with different value of α and the angle window of 15°.

The mean value as well as the standard deviation of the error across all patients is also shown in Table 2. This analysis is in line with our previous observations as controllers showed less error while incorporating the prediction element.

Table 2. The mean value as well as the standard deviation of the error across all patients.

	Binary	Adaptive PI	MRAC
Without prediction	−2.8 ± 0.94 mm	−2.7 ± 1 mm	−2.8 ± 0.89 mm
With prediction	−0.8 ± 1 mm	−0.5 ± 0.5 mm	−0.5 ± 0.7 mm

Remark: A homogeneous environment was considered for the simulations, primarily to enable accurate comparisons with the experiments performed in the tissue-mimicking phantom. We believe for a comparison study between different controllers a homogeneous environment can provide a valid comparison and does not affect the conclusions made in this paper. The effect of heterogeneities of tissue properties has been studied in our previous publications. The effect of changes in blood perfusion, ultrasound absorption coefficient and anatomical structures surrounding the prostate has been investigated in [10,11,18]. For example, we have observed that a reduced distance between the prostate boundary and surrounding bone can impact the amount of thermal dose accumulated in these structures [18]. High levels of blood perfusion can lead to an increase in undertreated volume fraction [18]. These changes have not affected the controllers’ stability; however, the performance can be degraded with variation in tissue parameters. For example as discussed by Goharrizi et al. [11] a low attenuation value or high perfusion value can cause some undertreatment for a large radius. We have also observed similar performance of the controller in pilot human studies in prostate which supports the robustness of this controller. One exception is the situation with large calcifications in the beam path, which obstruct penetration of the ultrasound beam and cause the controller to fail. Since this is a physical barrier, it will likely end up as an inclusion/exclusion criterion for these treatments. Where we believe future research efforts are required is in obtaining stable and accurate MR temperature measurements in the prostate gland and developing more advanced signal conditioning approaches to filter out any artefacts in the temperature measurement.

Comparison study

Three controllers were compared experimentally. Similar results were also obtained through simulations. Here, we only report the experimental results, since simulation results are consistent with experimental observations.

Experimental studies were conducted in a tissue-mimicking gel phantom (Zerdine, CIRS, Norfolk, VA). The ultrasound applicator was inserted into the phantom and attached to the MRI-compatible rotary motor. Imaging was performed using a 16-channel pelvic coil array on a 3-Tesla MRI (Achieva, Philips Healthcare, Best, Netherlands). The pixel-wise precision of the MR thermometry was approximately 0.5–1.0 °C in these experiments. We have also measured the precision of MR thermometry in the prostate gland in human volunteers and observed a similar uncertainty [19]. In our experience with experimental implementation of this technology, we can align the MRI slices with a precision of approximately one pixel, which is 1 mm. This corresponds to 20% of the element width, and results in acceptable performance. Philips Sonalleve® therapy planning software was used to run the experiments.

In these experiments the treatment was delivered over an angular sector that represented the posterior region of the prostate gland. This was done to mimic a targeted treatment of a localised tumour as was performed by our group in a parallel clinical study. The shape of the boundary was also selected to have sharp changes in radius to make the effect of heat diffusion more severe. The rate of change of radii in the clinical targets we have investigated was much slower, and within the envelope of acceptable performance for the controller.

The results of a 3D treatment considering two elements are shown in Figure 7. As can be seen, the radial targeting was very similar for all three controllers. A more quantitative view of the error between the targeted boundary and treated boundary has been plotted for each controller and each slice in Figure 8. The adaptive PI controller initially showed less error, but eventually all of the controllers turned out to have the same performance. We believe that the influence of the prediction element is very strong and therefore makes each controller perform equally.

Figure 7. The pattern of conformal thermal heating for two elements using different controllers (experimental results).

Figure 8. The difference between targeted and treated boundaries for two elements using different controllers.

In order to show the effectiveness of the prediction value we performed two more experiments with the binary controller using the extreme values of α. The results are shown in Figure 9. As one can see with α = 1, greater overtreatment was observed as predicted by the simulations. These results confirm the strong influence of the predictive element for this thermal treatment.

Figure 9. The pattern of conformal thermal heating for values of α using binary controller (experimental results).

Concluding remarks

In this paper a novel prediction methodology was developed to address the observation of overtreatment in MRI-controlled, minimally invasive ultrasound therapy. This prediction element could be integrated into any closed-loop control systems involving a moving ultrasound beam. The radial targeting accuracy of three controllers (adaptive PI, model reference adaptive, binary) was evaluated through simulation and experimental studies. The results were consistent and demonstrated equivalent performance between controllers. The dominant influence on radial targeting accuracy was the prediction element described in this paper. Based on these results, a binary controller with a predictive element may provide the best balance of performance and simplicity for this application.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.