Quantitative analysis of heat-source estimation of high-intensity focused ultrasound using thermal strain imaging

The temperature rise induced by HIFU irradiation results in temporal shifts of RF-echo signals. The temporal shift originates from two physical phenomena: thermal expansion of the propagation medium and local change in its speed of sound, which correspond to a physical and apparent shift, respectively. The temperature change $\delta \theta \left(z\right)$ is expressed with a linear coefficient of the thermal expansion $\alpha \left(z\right)$ and that of the temperature dependence of sound speed $\beta \left(z\right)\,\,$as

$$\begin{eqnarray}&&\delta \theta \left(z\right)=\displaystyle \frac{{c}_{0}\left(z\right)}{2}\left(\displaystyle \frac{1}{\alpha \left(z\right)-\beta \left(z\right)}\right)\cdot \displaystyle \frac{\partial }{\partial z}\left(\delta t\left(z\right)\right),\end{eqnarray} \tag{ 1 }$$

where ${c}_{0}\left(z\right)$ is the speed of sound at the initial temperature, $\delta t\left(z\right)$ is a time-shift on the echo at a depth z. The model can be applied to a limited temperature change within a range where the temperature dependence of sound speed is approximated to be linear. In addition, the coefficient $\alpha \left(z\right),\beta \left(z\right),$ and ${c}_{0}\left(z\right)\,$were assumed to be constant with respect to depth z in a uniform tissue. The typical value of $\beta \left(z\right)$ ranges from $0.7\times {10}^{-3}$ to $1.3\times {10}^{-3}\,/{\rm{K}}$ for water bearing tissue around body temperature, and is at least an order of magnitude larger than the coefficient of thermal expansion. ^20–22) All experiments used a tissue mimicking material (TMM) phantom produced based on the IEC60601-2-37, whose mass fractions and acoustic parameters are shown in Tables I and II. As for the TMM phantom, the coefficient $\beta \left(z\right)$ is assumed to be $1.3\times {10}^{-3}/{\rm{K}}$ in the temperature range of 23 °C–37 °C, calculated from the temperature coefficient of 2.1 m s⁻¹ K⁻¹ as reported on the acoustic properties of an IEC agar-based TMM phantom. ²³⁾ The coefficient $\alpha \left(z\right)$ is regarded as $0.1\times {10}^{-3}\,{{\rm{K}}}^{-1}$as an approximation of water ²⁴⁾ considering that 83% of the phantom is composed of water $\left[\alpha \left(z\right)=0.084\times {10}^{-3}\,{{\rm{K}}}^{-1}\right]$ ²⁵⁾ and 11% of glycerol [$\alpha \left(z\right)=0.16\times {10}^{-3}\,{{\rm{K}}}^{-1}$] ²⁵⁾ as shown in Table I.

Figure 2 shows the experimental setup for the measurement of the tissue-dependent coefficient $\alpha -\beta$ of the phantom. A rectangular (outer dimensions of 95 mm × 55 mm × 20 mm) phantom (OST TMM3S-125) was fixed in an acrylic tank filled with water. It was filled in a polyethylene film bag with zipper to prevent leaching of glycerol during a prolonged experiment. The temperature inside the phantom was varied at six different steps from approximately 20 °C–41 °C by heating the water in the tank with a throw type heater and stirring the heated water. At each temperature step, ultrasound imaging was performed when the temperature of the water measured by a thermometer (HORIBA OM-51) and that of inside the phantom measured by a thermocouple (CHINO SCHS1-0) almost coincided. The thermocouple was connected to a data logger (GRAPHTEC midi LOGGER GL900) and recorded at a sampling frequency of 500 Hz. For strain measurement, a phased array imaging probe (Hitachi Aloka Medical UST-52105, driven at a center frequency of 3.47 MHz) was placed 49 mm above the phantom surface and connected to a programmable ultrasound imaging system (Verasonics Vantage 256). RF data were acquired with 10 plane waves steered at angles ranged from −6° to 6° at a sampling frequency of 13.89 MHz. Ten consecutive in-phase/quadrature (IQ) data were spatially compounded in constructing as a frame. By applying a 2D autocorrelation method between the reference and each post-heating frame, five displacement maps were obtained from the phase difference at the center frequency determined based on the frequency spectrum of the IQ data. ^26,27) The parameter setting used for 2D autocorrelation is shown in Table III. The axial slope of the displacement was divided by the temperature rise measured by the thermocouple to calculate the tissue dependence coefficient $\alpha -\beta .$

Zoom In Zoom Out Reset image size

The temperature rise induced by HIFU irradiation was estimated using the measured values of $\alpha -\beta .$ As shown in the Fig. 3, experiments were performed in an acrylic tank containing degassed and deionized water (saturation concentration of dissolved oxygen of 20%–30%) at room temperature. A cylindrical (57 mm in diameter and 70 mm in height) phantom (OST TMM3S-0150D) was placed to contain the HIFU focal spot. The mass fractions and acoustic properties of the phantom are the same as those of the TMM3S-125 model described above. On the side wall of the tank, a 256-element two-dimensional-array HIFU transducer (120 mm in focal length and 120 mm in outer diameter) was installed, connected to a 128-channel staircase voltage driving system (Microsonic HIFU EU9144-2), and driven at a frequency of 1 MHz. In the central hole of the HIFU transducer, a phased array imaging probe (Hitachi Aloka Medical UST-52105) was fixed and controlled by a programmable ultrasound imaging system (Verasonics Vantage 256) to drive at a center frequency of 3.47 MHz and acquire RF echo signals.

Zoom In Zoom Out Reset image size

Figure 4 shows the sequence of HIFU exposure and RF data acquisition in this experiment.

Zoom In Zoom Out Reset image size

Prior to inducing the temperature rise, reference RF data were acquired with 5 plane waves steered at angles of −6°, −3°, 0°, 3°, and 6° at a sampling frequency of 13.89 MHz. After that, HIFU was irradiated to induce temperature rise of several degrees in Celsius, with a constant duration of 0.9 s at an intensity of 330 W cm⁻². To obtain frames to be compared, RF data were collected from 2 to 28 s after the end of HIFU irradiation. The consecutive IQ data at five angles were acquired and compounded to construct a frame. The displacement distributions were calculated in the same way as described in Sect. 2.2 with parameter setting shown in Table III. ^26,27) As shown in Fig. 5, each displacement map was spatially differentiated along the axial direction, and smoothed using a 1D Savitzky–Golay and a 2D Gaussian filter.

Zoom In Zoom Out Reset image size

In order to estimate the absorption coefficient of the phantom, the temperature rise induced by HIFU irradiation was directly measured with a thermocouple. The experiment was performed in the experimental setup described in Fig. 6. A thermocouple with a diameter of 0.3 mm was inserted at the point where the temperature rise with HIFU duration of 0.9 s was maximum. The temperature measured by the thermocouple was recorded at a sampling frequency of 500 Hz through the data logger. A low-pass filter with a cutoff frequency of 5 Hz was applied to reduce the noise. The method used for data analysis is explained in the next subsection.

Zoom In Zoom Out Reset image size

With regard to the measurement of temperature rise shown in Fig. 6, a tip of thermocouple was placed in an ultrasound field, so the artifact associated with viscous heating was no longer negligible, which occurred at an interface between a thermocouple and the surrounding tissue by the relative motion between them. ²⁸⁾ The time constant of thermal diffusion of viscous heating is known to be much shorter than that of heating due to ultrasonic absorption of tissue, and the contribution of the latter to the temperature rise becomes dominant certain time after the HIFU irradiation. In order to separate the temperature rise into the ultrasonic absorption and viscous heating components, the measured temperature rises were fitted by using the least-squares method. ^29,30) The fitting model of temperature rises were simulated by HITU simulator, where a wide-angle Khokhlov–Zabolotkaya–Kuznetsov and bioheat transfer equations are solved for the pressure and temperature distribution, respectively. ³¹⁾ The computational domain was divided into three layers along the depth, and the phantom (Material 2) was assumed to be sandwiched between two water domains (Materials 1 and 3) with parameter settings as shown in Table IV. The fraction of attenuation due to absorption in water was set to be 0 to neglect the heat generation in water. For a viscous heating model, a minute heat source of cylinder with 0.2 mm in diameter and 0.2 mm in height was assumed at the position of the thermocouple. The percentage of absorption in attenuation in the TMM phantom, which is shown in material 2 in Table IV was varied so that the standard deviation in four results with HIFU durations of 0.2, 0.4, 0.6, and 0.9 s should be the smallest.

Figures 7(a)–7(e) show displacement distribution maps with five different temperature conditions. The temperature rise $\delta \theta$ = 3.9 °C, 9.2 °C, 12.3 °C, 17.5 °C, and 21.5 °C from left to right. In the displacement maps, the imaging ultrasound was irradiated downward, and positive displacement indicates the upward shift. The slope of the displacement was calculated for each line in the depth direction for a portion of the depth as shown in Figs. 7(f)–7(j), and averaged value was taken as the thermal strain at the temperature rise. Figure 8 shows the measured thermal strain as a function of temperature rise inside the phantom for three different samples. A red line represents a linear approximation passing through zero. According to the slope of the fitting line, the averaged tissue-dependent coefficient was estimated as $\alpha -\beta =\,-1.38\,\times {10}^{-3}\,{{\rm{K}}}^{-{\rm{1}}}$ within the range of measurement, which was close to the literature values of $\alpha =\,0.1\times {10}^{-3}\,{{\rm{K}}}^{-{\rm{1}}}$(water) and $\beta =1.3\times {10}^{-3}$ K⁻¹ (TMM phantom).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the temperature rise measured by the thermocouple with HIFU durations of 0.2, 0.4, 0.6, and 0.9 s (n = 3). The red dash and black solid lines represent the measured temperature rise and the fitting result of the simulation, respectively. The green and blue solid lines show the temperature rise components from ultrasonic absorption and viscous heating, respectively, the sum of which gives the fitting result. From the result that the fraction of absorption in attenuation was 0.9 with the smallest standard deviation of the coefficient, the absorption coefficient of the TMM phantom was estimated to be 0.45 dB cm⁻¹ MHz⁻¹ in this experiment. The standard deviation of the coefficient was 7.9% for four HIFU durations. The subsequent simulations of temperature rise were conducted using the absorption coefficient of 0.45 dB cm⁻¹ MHz⁻¹.

Zoom In Zoom Out Reset image size

Figure 10(a) shows a set of displacement distribution maps acquired at 2–28 s after the end of HIFU irradiation, and Fig. 10(b) shows corresponding cross-correlation maps. The displacement at and below the focal point has positive value, which means that the tissue in the region was displaced upward. As shown in Fig. 10(b), the correlation was lower at and behind the heated region. Although the correlation is the lowest at 2 s from the end of HIFU irradiation, it is improved to more than 0.995 at 28 s from the end of HIFU irradiation. Figure 10(c) shows a set of thermal strain maps acquired at 2–28 s after the end of HIFU irradiation. The position with a maximum absolute value of thermal strain agreed well with the geometric focal depth of 70 mm in the images. Moreover, it is observed that the thermal strain distribution slightly expanded to the front and lateral regions with lapse of time due to the thermal conduction.

Zoom In Zoom Out Reset image size

In order to quantitatively evaluate the temperature rise induced by HIFU irradiation, the temperature rise was obtained from the measured thermal strain using the measured value of $\alpha -\beta$ into Eq. (1). Figure 11 shows a comparison between the estimated and simulated temperature rises over time from the end of HIFU irradiation. The estimated temperature rise (n = 9) was averaged in ROI (Lateral: −3.4–6.4 mm, Axial: 65–70 mm), and simulated temperature rise was also averaged in the corresponding ROI. The black solid line represents the estimated temperature rise with measured value of $\alpha -\beta .$ The error bars show the standard deviation. As for simulation results, the red solid line shows the temperature rise averaged in the same ROI for the estimation, and the blue solid line shows the temperature rise considering the point spread function (PSF) in the elevational direction. As shown in Fig. 12, the simulated temperature was convoluted with the PSF in the elevational width of 3.5 mm, which is close to the full width half maximum of the elevational PSF. The estimated temperature rise based on the measured value of $\alpha -\beta$ shows good agreement with the simulated temperature rise. The agreement was even better when the spatially averaged effect owing to the size of PSF was considered.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size