High intensity focused ultrasound (HIFU) is becoming a well-accepted ablation therapy because of its patient quick recovery, focal treatment and non-invasiveness.1 To ablate a target tissue, thermal monitoring is required to identify the ablation completeness. Moreover, thermal monitoring helps to prevent destroying surrounding healthy tissues. Many medical imaging modalities have been studied to monitor the temperature during ablation therapy such as MRI,2 periodic X-rays3 and medical ultrasound.4, 5 Medical ultrasound (US) imaging is a cost-effective technique and its accessibility is higher than the other medical imaging modalities. For these reasons, various ultrasound thermometry methods have been proposed.
Temporal ultrasound tomography method can detect the speed of sound (SOS) change caused by the temperature increase. Ultrasound travels through tissues at different velocities when the temperature is not uniform in the propagating path. This phenomenon enables ultrasound thermal monitoring during an ablation procedure.6 The different SOS in a region of interest can be recovered by solving a minimization problem using recorded time-of-flight (TOF) information7 in order to track the temperature changes.
Thermal HIFU ablation systems utilize focused ultrasound waves to destroy target tissues by increasing their temperature.8 By adding a receiving element, the SOS in the ablation area can be estimated using TOFs from the HIFU elements. We showed the feasibility of the HIFU-US thermal monitoring method using an external ultrasound element in previous work.9
In this paper, we compare two methods to reconstruct SOS volume: an iterative optimization and a batch optimization method. In addition, the location of the ultrasound element affects the estimation accuracy. Indeed, the paths between each HIFU element and the receiving ultrasound element bring valuable information if they intersect with the ablation area. Their intersection lengths with the ablation area depend on the location of the ultrasound element. For this reason, we analyze the effect of the receiving ultrasound element location. Finally, as the accuracy is also affected by the noise level, we analyze if this error can be minimized by choosing an appropriate receiving ultrasound element position.
HIFU-US THERMAL MONITORING
We add an ultrasound element to collect ultrasound pressure waves transmitted by each HIFU element. The TOF information is used to recover the SOS in a volume of interest. As the SOS is related to the tissue temperature, this method can be used for thermal monitoring purpose.
Ultrasound thermal monitoring setup
Our ultrasound thermal monitoring method requires an external ultrasound receiving element, which can be integrated with the HIFU ablation system without the need of any additional ultrasound probes since the HIFU system is designed to transmit ultrasound. We have fabricated a small size ultrasound element able to receive the ultrasound pulses from a wide angle.
Figure 1 shows an overview of our HIFU ultrasound thermal monitoring setup. The receiving element is placed on top of the tissue.
A MR-HIFU system (Sonalleve V2, Profound Medical, Toronto, Canada) is simulated for this study. The HIFU system consists of 256 elements, and their arrangement is illustrated in Figure 1. We define the (x,y,z) coordinate of HIFU coordinate system in mm unit. Natural focus of the HIFU system is the origin of the coordinate system. The direction orthogonal to the patient bed is defined as x axis, and y and z axis are defined as the transverse and longitudinal direction of the bed.
In our setup, the HIFU-US thermal monitoring system interleaves ablation and monitoring phases. Continuous waves focus on the target during an ablation phase, and then each HIFU element transmit short ultrasound pulses one after each others. The TOFs from each HIFU element are acquired sequentially in the monitoring phase.
A localization of the element is required to order to relate the measured TOFs with SOS changes. It is performed before ablation, during the first monitoring phase. The 3D coordinates of each HIFU element are known, so we can localize the ultrasound element coordinate with respect to the HIFU coordinate system using a minimum of three TOF datasets which are not non-linearly aligned in 3D space. We recover the SOS in a region of interest (ROI) using TOFs change in the monitoring phase.
Speed of sound reconstruction
We estimate the speed of sound in each voxel in the volume of interest. As we collect the TOFs between every HIFU element and the ultrasound receiving element, we have as many equations as the product of the number of HIFU elements by the number of the ultrasound receiving element. This number of equations is smaller than the number of unknown SOS in each voxel. Therefore, we use a computational model to have an priori knowledge of the temperatures. The computational model is described in the following section 2.3. We define as a layer a group of voxels which are expected to have the same temperature according to the HIFU simulation. Thus we are able to estimate the SOS in each layer with the limited number of equations available. Finally, the SOS is converted to 3D temperature images.
Equation 1 presents the optimization problem.9 Unknown x is a vector that contains the inversed SOS in each voxel. ToF vector is created from the segmented TOF between the HIFU elements and the ultrasound element. The system matrix P contains the intersecting length values in each voxel, which has a size of the number of HIFU by the number of voxels. The equality and inequality conditions are derived from computational model of HIFU.
Two optimization methods are employed to reconstruct a SOS image volume. First, a batch optimization process, where we minimize the objective function to solve for the SOS in the entire volume of interest simultaneously. Then, we implemented an iterative method to recover the SOS starting from the outer layers which usually intersect with more paths. Each acoustic propagation path intersects with a layer or multiple layers in different lengths. In the iterative optimization mode, we solve for the SOS starting from the outer layers considering only the TOF sets having the largest fractional lengths through it. We choose N largest sets to solve for the SOS in each layer. If the number of equation sets which are intersecting with the target layer is less than N, then we choose all sets which have non-zero intersecting value with this layer. Moreover, depending on the ultrasound element location, there could be some layers which does not intersect with any propagation path. In this case, the SOS in such layers are solved by the inequality constraint with the neighboring layers.
Computational modeling of HIFU
The physics-based HIFU simulation model includes both acoustic pressure simulation and heat propagation. It is based on nonlinear ultrasound propagation using a k-space model coupled with heat propagation in biological tissue using a reaction-diffusion equation, implemented using the Lattice Boltzmann Method (LBM). The acoustic simulation is performed first in order to simulate the ultrasound pressure and velocity fields, which are used as input of the heat propagation simulation. 3D longitudinal thermal volumes images are simulated every 1s of the procedure.
SIMULATION AND RESULT
Temperature and speed of sound images
One HIFU ablation consisting of three consecutive phases as described in more detail in another paper10 was used. One phase consists of a heating period of 20 seconds followed by a cooling period of 24 seconds. During the heating phase, all the HIFU elements transmit continuous waves at 1.2 MHz with a power of 50 W. The 256 HIFU elements are placed asymmetrically11 in the coordinates between (-140,-80,-80) and (-120,80,80), the focal heating point defining the origin. The temperature images generated by the HIFU computational model during the third and last cooling phase are averaged in one single image shown in figure 2. This averaged temperature image was used to generate a SOS volume using a converting equation specific for 2% agar and 2% silicon-dioxide material.9 This SOS image is used as ground truth in all the experiments described in this paper. The temperature ranges from the background temperature: 21.0°C, which is the room temperature to 36.4°C. It corresponds to a SOS range of 1500m/s to 1539m/s.
We grouped voxels in 25 different layers, defined with a step size of 0.6°C, thus every voxel with a the temperature below 21.6° C is considered as a background voxel. In other words, the ROI consists of the voxels with a temperature above 21.6°C. And to reduce the computation time, only the voxels in the ROI are considered. As the image resolution was 1mm isotropic, the number of voxels in the ROI was 18,121.
The ultrasound element location
The distribution of the intersection lengths between the ablation area and each of the path between one HIFU element and the receiving element changes with the position of the receiving ultrasound element. This distribution directly affects the SOS reconstruction.
In this simulation study, the ablation is performed 25 mm beneath the surface. We assume that the tissue surface is flat and so the element can be located anywhere on the surface plane. We reconstructed the SOS volume image for different receiving element locations between (25,-10,-10) and (25,10,10) at 2 mm step in both y and z axis. In each case, both the batch optimization and the iterative optimization are compared. For the iterative approach, different numbers of N (2, 8, 32, 128, and 256) are used. Results are shown in Figure 3. The estimation error is inherent to the approach since we group the voxels using a temperature step of 0.6°C. As we increase N, the error generally decreases since the method utilizes more information from more TOF sets. However, lower N iterative approach performs better at certain locations. It depends on the intersection length distribution profile. For both the iterative and batch optimizations, the area between (25,-2,-2) and (25,2,2) has largest errors since all the ultrasound pulse propagation paths do not go through the heating center (at least the 5 more central layers) as shown in the Figure 4.
The intersection length distributions for different ultrasound element locations are shown in Figure 4. The intersection length distribution with a triangular profile and including the inner layers yield less error while the distributions with a more flat profile and not including the inner layers result in larger estimation error. Although most of the propagating paths pass through the ablation area for each element location, if the TOF affected by the inner layers is not acquired then it results in relatively larger estimation errors. At location b: (25, -10, 8), almost half of the TOF does not intersect with the ablation area. However, the acoustic paths intersect with every layer of the ROI and therefore the SOS reconstruction is more accurate.
All the computations were performed with a work station with i7-4790 CPU @ 3.6GHz, and 16GB Ram. A Matlab linear least square solver is used to minimize Equation 1. In the iterative approach, the computation could be done in around 1 second in the most of cases when N is less than 8, and the computation time increased with N. When N reaches 32, the computation time becomes similar to the one of the batch optimization. On average for the 121 different element locations considered, the computation time of the batch optimization was 2.1 seconds. The iterative approach with lower N could suffer from relatively larger errors, but the execution time is faster than for the batch approach.
We reconstructed SOS images using two different methods: the batch optimization and the iterative approach. To compare the sensitivity of the two methods to the presence of noise in the acquired datasets, we fix the receiving element to a location providing similar errors for both methods. When the ultrasound element is located at (25,-4,-8), the batch optimization maximum SOS error was 1.03m/s, and the iterative method had an error of 1.09m/s. The mean error in ROI were 0.35m/s and 0.36m/s respectively.
We induced a Gaussian random noise within ±5us on the TOF datasets. This noise corresponds to a segmentation error of up to 5 cycles of the ultrasound pulse when we use a 1 Mhz transmitting pulse. 100 simulations were performed with this random noise. Results are shown in Figure 5. The standard deviation of the maximum error with the batch optimization was 0.230m/s while it was 0.208m/s using the iterative approach (N=32). The median error was less than 1.5m/s for both methods, which corresponds to less than 1°C error in temperature at 30° C. With the Gaussian noise we used in this study, the maximum error range was more stable with the iterative method.
DISCUSSION AND CONCLUSION
We introduce two speed of sound reconstruction methods for the temperature monitoring during HIFU ablation and study the effects of both the receiving ultrasound element location and noise in the TOF datasets. Equipment such as the HIFU system and the sampling devices can induce system delays. Moreover, electronic and acoustic noises can also affect the accuracy of the TOF segmentation.
The least square method is generally sensitive to noise and outliers. A noise in the TOF dataset can considerably affect the estimation accuracy if the flight-of-time of an intersection length is relatively shorter compared to the noise level. When a noise is induced in a path on a certain layer, the error could also affect the other layers. Especially, a wrong estimate in an outer layer can affect the inner layer SOS solution to a large extent.
The batch optimization method and the iterative SOS reconstruction method result in different errors. The simulation results show that the iterative method is feasible for the SOS reconstruction. When we are dealing with more noise in the TOF data, the iterative method may be preferred since it often shows less error deviation.
Clinical applications have to deal with patient motions, which would be a crucial factor of the monitoring accuracy. In our settings, when the system detects sudden change in TOFs, the thermal monitoring algorithm can restart from the localization phase and recover the temperature from prior SOS information.
The proposed SOS reconstruction approach can be further improved. First, the error would be decreased by increasing the number of ultrasound elements so that more TOFs from different locations could be acquired simultaneously. As the number of TOF data increase, we can divide the layers in finer steps which could increase the estimation accuracy, which is limited by the layer definition based on a heat propagation simulation model. The computation time can also be further decreased with a high-performance workstation. In the current implementation, the optimization method generated the SOS volume in 1 second.
Ultrasound thermal monitoring has the advantages of its portability and accessibility. The proposed SOS reconstruction methods can be used for HIFU-US thermal monitoring by adding minimal devices using batch and iterative optimization methods. We studied the noise sensitivity and the effect of the ultrasound element location on the estimation accuracy. Using the simulation results, we can optimize the proposed HIFU ultrasound thermal monitoring method.
Research reported in this paper was supported by the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health under award number 01EB021396. and National Science Foundation under proposal number 1653322.
Brace, C. L., Mistretta, C. A., Hinshaw, J. L., and Lee, F. T., “Periodic contrast-enhanced computed tomography for thermal ablation monitoring: A feasibility study,” Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 4299–4302 (Sept 2009).Google Scholar
Norton, S. J., Testardi, L. R., and Wadley, H. N. G., “Reconstructing internal temperature distributions from ultrasonic time-of-flight tomography and dimensional resonance measurements,” Ultrasonics Symposium, 850–855 (Oct 1983).Google Scholar
Kim, Y., Guo, X., Cheng, A., and Boctor, E. M., “Speed of sound estimation with active pzt element for thermal monitoring during ablation therapy: feasibility study,” Proc. of SPIE Medical Imaging 9790, 97901K–97901K-8 (2016).Google Scholar
Kim, Y., Audigier, C., Ellens, N., and Boctor, E. M., “A novel 3D ultrasound thermometry method for HIFU ablation using an ultrasound element,” IEEE International Ultrasonics Symposium (IUS), 1–4 (Sept 2017).Google Scholar
Audiger, C., Kim, Y., Ellens, N., and Boctor, E. M., “Simulation of high intensity focused ultrasound ablation to enable ultrasound thermal monitoring,” Proc. SPIE Medical Imaging (2017).Google Scholar
Mougenot, C., “L’ asservissement par IRM dun reséau matriciel ultrasonore et ses applications therapeutiques,” PhD thesis, L’Université de Bordeaux (2005).Google Scholar