Physics-Based Simulation to Enable Ultrasound Monitoring of HIFU Ablation: An MRI Validation

Abstract

High intensity focused ultrasound (HIFU) is used to ablate pathological tissue non-invasively, but reliable and real-time thermal monitoring is crucial to ensure a safe and effective procedure. It can be provided by MRI, which is an expensive and cumbersome modality.

We propose a monitoring method that enables real-time assessment of temperature distribution by combining intra-operative ultrasound (US) with physics-based simulation. During the ablation, changes in acoustic properties due to rising temperature are monitored using an external US sensor. A physics-based HIFU simulation model is then used to generate 3D temperature maps at high temporal and spatial resolutions. Our method leverages current HIFU systems with external low-cost and MR-compatible US sensors, thus allowing its validation against MR thermometry, the gold-standard clinical temperature monitoring method.

We demonstrated in silico the method feasibility, performed sensitivity analysis and showed experimentally its applicability on phantom data using a clinical HIFU system. Promising results were obtained: a mean temperature error smaller than $1.5\,^\circ {\text {C}}$ was found in four experiments.

1 Introduction

In the past two decades, high intensity focused ultrasound (HIFU) has been used with generally good success for the non-invasive ablation of tumors in the prostate, uterus, bone, and breasts [1], along with the ablation of small volumes of neurological tissue for the treatment of essential tremor and Parkinson disease [2]. Even though its use has considerably expanded, a major limitation is still the lack of detailed and accurate real-time thermal information, needed to detect the boundary between ablated and non-ablated zones. The clinical end-point being to ensure a complete ablation while preserving as much healthy tissue as possible. This kind of information can be provided with $\pm 1\,^\circ {\text {C}}$ accuracy by MRI [3], routinely used to guide HIFU in clinical settings [4].

However high temporal and spatial resolutions are needed to accommodate with the small and non-uniform ablation shape and to detect unexpected off-target heating. This is challenging to achieve over a large field of view with MRI as the scan time limits the volume coverage and spatial resolution. Typically, one to four adjacent slices perpendicular and one slice parallel to the acoustic beam path, with a voxel size of 2$\,\times \,$2$\,\times \,$5 mm, at 0.1–1 Hz are acquired. Moreover, this modality is expensive, cumbersome, and subject to patient contraindications due to claustrophobia, non-MRI safe implants or MR contrast material.

Compared to MRI, ultrasound (US) offers higher temporal and spatial resolutions, low cost, safety, mobility and ease of use. Several heat-induced echo-strain approaches relying on successive correlation of RF images, have been proposed for US thermometry [5]. Despite good results in benchtop experiments, they suffer from low SNR, uncertainties in the US speckles, weak temperature sensitivity and have failed to translate to clinical applications, mainly due to their high sensitivity to motion artifacts, against which the proposed method is robust since it relies on direct measurements in the microsecond range.

We present an inexpensive yet comprehensive method for ablation monitoring, that enables real-time assessment of temperature and therefore thermal dose distributions, via an integrative approach. Intra-operative time-of-flight (TOF) US measurements and patient-specific biophysical simulation are combined for mutual benefits. Each source of information alone has disadvantages. Accurate simulation of an ablation procedure requires the knowledge of patient-specific parameters, which might not be easy to acquire [6]. US thermometry alone is not robust enough to fully meet the clinical requirements for assessing the progression of in-vivo tissue ablation. We propose to leverage conventional HIFU system with external low-cost and MR-compatible US sensors to provide in addition to ablation, real-time US temperature monitoring. Indeed, as HIFU deposits acoustic thermal dosage, invaluable intra-operative information is usually omitted. With rising temperature and ablation progression, acoustic properties such as the speed of sound (SOS) and attenuation coefficient vary. This affects the TOF carried by the US pressure waves going through the ablation zone and propagating to the opposite end, which we intend to record by simply integrating US sensors. Moreover, the proposed approach allows to use MRI for validation.

US thermometry through tomographic SOS reconstruction from direct TOF measurements has previously been proposed [7]. But for HIFU monitoring, the tomographic problem is ill-posed. First, it is rank deficient as the acquired TOF data is sparse: equal to the number of HIFU elements: 256 with common clinical HIFU system, times the number of sensors employed. However, we aim to reconstruct the temperature at the voxel level. Moreover, the relationship between SOS and temperature is tissue-specific and linear only until a certain point (around 55–60 $^\circ {\text {C}}$). To tackle these, we propose to incorporate prior knowledge of biological and physical phenomena in thermal ablation through patient-specific computational modeling. 3D thermal maps at a high temporal and spatial resolution as well as TOF variations during ablation progression are simulated.

In this paper, we introduced the proposed method, presented simulation experiments and sensitivity analysis to evaluate its feasibility. In vitro validation against MRI in 4 phantom experiments were also performed.

2 Methods

HIFU ablation consists in the transmission of high intensity US pressure waves by all the HIFU elements. Each wave passes through the tissue with little effect and propagates to the opposite end. At the focal point where the beams converge, the energy reaches a useful thermal effect. By integrating external MR-compatible US sensors placed on the distal surface of the body (Fig. 1A), our method records invaluable direct time-of-flight (TOF) information, related to local temperature changes. Therefore, a large US thermometric cone defined by the HIFU aperture and the US external sensor is covered (Fig. 1B).

Fig. 1.

The system setup (A) a 3D rendering showing the HIFU system embedded in the patient bed, external MR-compatible ultrasound (US) sensors, and MRI gantry. (B) a zoomed-in schematic diagram highlights individually controlled HIFU elements, US sensor, both HIFU transmit cone and US thermometry cone, and the sagittal MR imaging plane. (C) timing diagram shows both HIFU ablation and monitoring phases.

2.1 Biophysical Modeling of HIFU Thermal Ablation

HIFU is modeled in two steps, ultrasound propagation followed by heat transfer.

Ultrasound Propagation: The nonlinear US wave propagation in heterogeneous medium is simulated based on a pseudo-spectral computation of the wave equation in k-space [8]. The US pressure field $p(\mathbf {x},t)$ [Pa] is computed during the ablation and monitoring phases.

Heat Source Term: $Q(\mathbf {x},t)$ [W/$m^3$] is computed from the US pressure as [9]:

$$\begin{aligned} Q(\mathbf {x},t) = \frac{\alpha f}{\rho _t c}|p(\mathbf {x},t)|^2 \end{aligned}$$

(1)

where $\alpha$ [Np/(m MHz)] is the acoustic absorption coefficient, f [MHz] the HIFU frequency, $\rho _t$ [kg/m$^3$], the tissue density and c [m/s] the speed of sound.

Heat Transfer Model: From a 3D anatomical image acquired before ablation, the temperature $T(\mathbf {x},t)$ evolution is computed in each voxel by solving the bioheat equation, as proposed in the Pennes model [10]:

$$\begin{aligned} \rho _t c_t \frac{\partial T(\mathbf {x},t)}{ \partial t} = Q(\mathbf {x},t) + \nabla \cdot ( d_t T(\mathbf {x},t)) + R(T_{b0}-T(\mathbf {x},t)) \end{aligned}$$

(2)

where $c_t$ [J/(kg K)], $d_t$ [W/(m K)] are the tissue heat capacity and conductivity. $T_{b0}$ is the blood temperature. R is the reaction coefficient, modeling the blood perfusion, which is set to zero in this study, as we deal with a phantom.

2.2 Ultrasound Thermal Monitoring

From the Forward Model: The biophysical model generates longitudinal 3D temperature maps, which can be used to plan and/or monitor the ablation. They are also converted into heterogeneous SOS maps given a temperature-to-SOS curve (Fig. 2 shows such a curve measured in a phantom). Thus, US wavefronts are simulated as emitted from each HIFU element and received by the US sensor with temperature-induced changes in their TOF. For monitoring, the recorded TOFs are compared to those predicted by the forward model: if they are similar, then the ablation is going as expected and the simulated temperature maps can be used. If the values diverge, the ablation should be stopped. Thus, insufficient ablation in the target region and unexpected off-target heating can now be detected during the procedure.

From tomographic SOS Reconstruction: During each monitoring phase $t_{m}$, 3D SOS volumes are reconstructed by optimizing Eq. 3 using the acquired TOF, provided that the number of equations is at least equal to the number of unknowns, i.e. the SOS in each voxel of the 3D volume. As the acquired TOF data is limited, additional constraints are needed. To reduce the number of unknowns, we created layer maps $M_{t_{m}}$. Each $M_{t_m}$ includes $N_{t_m}$ different layers grouping voxels expected to have the same temperature according to the forward model.

$$\begin{aligned} \begin{aligned} \underset{x}{\text {min}} \Vert Sx - TOF_{aquired} \Vert ^2 \quad \text {subject to } \quad&A_{eq}\cdot x = b_{eq}, \\&A_{ineq}\cdot x \le b_{ineq}, \\&SOS_{min} \le 1/x \le SOS_{max} \\ \end{aligned} \end{aligned}$$

(3)

where the vector x represents the inverse of the SOS, the matrix S contains the intersection lengths through each voxel for the paths between the HIFU elements and the US sensor, $\mathcal {P_{HIFU \rightarrow US}}$. Constraints between voxels in the same layer ($A_{eq}$, $b_{eq}$) and different layers ($A_{ineq}$, $b_{ineq}$) are computed based on $M_{t_m}$. The solution is also bounded by a feasible SOS range. From the estimated SOS, the temperature can be recovered using a temperature-to-SOS curve (Fig. 2).

3 Experiments and Results

We used a clinical MR-HIFU system (Sonalleve V2, Profound Medical, Toronto, Canada) providing in real-time 2 MR temperature images in the coronal and sagittal planes for comparison. We reprogrammed the HIFU system [11] to perform three consecutive cycles consisting of a heating phase at 50 W and 1.2 MHz (all elements continuous wave) for 10 s and a monitoring phase with an element-by-element acoustic interrogation at 2 W (40 cycle pulses, 128 elements sequentially) for 24 s to determine TOF (Fig. 1C). We fabricated a receiving US sensor, made of a 2.5-mm diameter tube of Lead Zirconate Titate material (PZT-5H). A 3-m long wire is connected to an oscilloscope located outside of the MR room. The US sensor was placed on top of the phantom, about 15 cm from the transducer and was localized using TOF measured prior to ablation. The sensor produced negligible artifacts on the MR images (Fig. 2, left).

Fig. 2.

(Left) experimental setup: HIFU is performed on a phantom under MR thermometry. The MR-compatible US sensor is placed on top of the phantom to acquire time-of-flight (TOF) data during the monitoring phase of the modified protocol. (Right) the phantom-specific temperature-to-SOS curve.

3.1 Sensitivity Analysis

The protocol described above was simulated using the forward model with both the ablation and monitoring phases. The US propagation and heat transfer are computed on the same Cartesian grid including a perfect match layer (PML = 8), with a spatial resolution of $\delta x= 1.3$ mm, different time steps: $\delta t= 12.5\,{\upmu }$s and $\delta t = 0.1$ s respectively, and using the optimized CPU version of k-Wave 1.2.1¹. Temperature images of 1.3$\,\times \,$1.3$\,\times \,$1.3 mm were generated every 1 s by the forward model. Layer maps $M_{t_m}$ were generated with a temperature step of 0.6 $^\circ {\text {C}}$. For example, at the end of the third ablation phase $t_3$, when the maximal temperature is reached, a map of $N_{t_3}$ = 26 layers was generated.

Effect of US Element Location: As the US sensor defines the US thermometry cone, its location with respect to the HIFU system and to the heated region highly affects the monitoring temperature accuracy. To study this effect, we simulated 4 different US sensor locations at the phantom top surface, 1 cm away from each other in each direction (the same setting was replicated in the phantom experiment, as detailed below). For each sensor location, we computed a matrix I, made of the cumulative length of the intersection of $\mathcal {P_{HIFU \rightarrow US}}$ with each of the layers of $M_{t_3}$. From the example illustrated in Fig. 3, it can be observed that at location A and D, most of the layers are covered by several paths. However, at location B, most of the paths from the central HIFU elements are not going through any layers, making it more difficult to reconstruct accurately SOS maps. This is also true for location C although to a lesser extent.

Fig. 3.

(Left) one path from one HIFU element to the US receiver going through the layer map is displayed over a simulated temperature image. 4 layers with $a< b< c < d$ are shown. (Right) matrices I made of the cumulative length of the intersection of $\mathcal {P_{HIFU \rightarrow US}}$ with each of the layers in $M_{t_3}$ are shown for 4 US sensor positions.

Fig. 4.

(Left) MR (top) and simulated (bottom) thermal images compare qualitatively well in a ROI of 75$\,\times \,$75 mm, centered around the targeted region. (Right) measured (top) and simulated (bottom) TOF from HIFU element 241 to the US sensor in position A. Delays of 0.1, 0.2, 0.3 ${\upmu }$s occur after the first, second and third heating phases.

Effect of the Number of US Elements: Multiple US sensors receiving simultaneously TOF at different locations can be used to improve the method accuracy, as the matrix I is less sparse. As it can be observed in Fig. 3, a better sampling of the layers (less zeros in I) can be achieved with certain combination of the 4 sensors. For example, by combining US sensors at location A and B, information about the outer layers 23 to 25, is obtained with high cumulative lengths by the HIFU elements 36 to 112 from the US sensor at location A, whereas information about the inner layer 1 comes from location B.

Fig. 5.

Error between the temperature estimated by the SOS reconstruction algorithm and the coronal MR image. The TOF from 1, 2, 3 or 4 sensors are used. The mean error in yellow is lower than 1 $^\circ {\text {C}}$ in each case. The max error at $t_1$, $t_2$ and $t_3$ appears in blue, cyan and pink. The overall max error decreases when we increase the number of elements, as shown by the black horizontal lines.

3.2 Phantom Feasibility Study

Four experiments with different US sensor locations were performed on an isotropic and homogeneous phantom made of 2%-agar and 2%-silicon-dioxide. Its specific temperature-to-SOS curve was measured pre-operatively (Fig. 2). First, the acquired MR thermal images and the ones generated by the forward model were compared. As shown in Fig. 4 at $t_3$, temperature differences in a ROI of 75$\,\times \,$75mm were 0.7 ± 1.2 $^\circ {\text {C}}$ and 1.6 ± 1.9 $^\circ {\text {C}}$ on average, with a maximum of 6.7 $^\circ {\text {C}}$ and 11.7 $^\circ {\text {C}}$ in the coronal and sagittal planes, respectively. TOF simulated at baseline and after the first, second and third heating phases were in agreement with the measures (Fig. 4). The delays caused by the temperature changes were computed by cross-correlation between signals received before and during ablation.

To evaluate in vitro the effect of multiple sensors acquiring TOF simultaneously, individual measurements obtained sequentially at the 4 different locations were grouped to mimic the monitoring by 2, 3 or 4 sensors. It was possible since we waited for the phantom to return to room temperature between each experiment. We analyzed all the combinations using 1, 2, 3 and 4 sensors. In each of the 15 scenarios, temperature images are generated using the tomographic SOS reconstruction and compared to MRI. As illustrated in Fig. 5 for the coronal plane, the algorithm accuracy highly depends on the position and number of the US sensors, as predicted in the above sensitivity analysis study. The overall max error decreases with the number of US elements employed. Similar results were obtained in the sagittal plane, with a mean error lower than 1.5 $^\circ {\text {C}}$ in each case.

4 Discussion and Conclusion

In this work, we used nominal parameters from the literature, but the biophysical model handles the presence of blood, different tissue types, and can be personalized to simulate patient-specific ablation responses to improve its accuracy in vivo [6]. As the simulation runs fast, model parameters could be personalized from intra-operative measurements during a first ablation phase and then used in the following phases. Different temperature-to-SOS curves could also be used.

As the heating is paused during monitoring, this period is desired to be as short as possible. To be more effective, one could minimize the switching time while cycling through the HIFU elements without inducing cross-talk. One could also sonicate pulses from multiple elements at once and deconvolve the received signals geometrically. Finally, we could investigate whether a smaller subset of the HIFU elements could be sufficient for temperature monitoring.

In conclusion, we have shown that biophysical model simulating the effect of treatment on patient-specific data, can be combined with US information directly recorded from HIFU signals to reconstruct intra-operative 3D thermal maps. This method demonstrated low temperature error when compared to MRI. While this work is a proof of concept with simulation and preliminary but solid phantom results, in vivo experiments are warranted to determine the viability of this US thermal monitoring approach. It promises to increase the safety, efficacy and cost-effectiveness of non-invasive thermal ablation. By offering an affordable alternative to MRI, it will for example transform the treatment of uterine fibroid into an outpatient procedure, improving the workflow of gynecologists who typically diagnose the disease but cannot perform MR-guided HIFU. By shifting the guidance to US, this procedure will be more widely adopted and employed.