Modelling the temperature evolution of bone under high intensity focused ultrasound

The model takes the actual geometry of the phased array transducer used in the Philips Sonalleve system into account (Köhler et al 2009). This phased-array HIFU transducer consists of 256 elements placed on a spherical shell (12 cm radius of curvature, 13 cm aperture) operating at 1.4 MHz, creating a cigar-shaped focus point with dimensions of approximately 2 mm width and 7 mm length. The flow of acoustic power from each transducer element is described by a number of rays, that have as their target randomly selected points within a 2 mm diameter circle in the focal plane, approximating the focal region described above. Each ray has an initial power ${{P}_{\text{tot}}}/{{N}_{\text{rays}}}$, where ${{P}_{\text{tot}}}$ is the total acoustic power and ${{N}_{\text{rays}}}$ (>=100 000) is the total number of rays starting from all the transducer elements. In one material type the rays propagate along straight lines and their power is attenuated using the attenuation coefficient of the considered material. If I(s) is the power flux in a ray after s cm in a material with the attenuation coefficient (for the pressure) α, then, since the acoustic power is proportional to the square of the pressure, $I(s)=I(0){{\text{e}}^{-2\alpha s}}$.

Waves propagate differently in liquids and in solids. In liquids only longitudinal (pressure) waves can travel, whereas in solids shear waves are also possible. For our model, muscle and marrow are considered as being liquid, with ultrasound properties similar to water. Thus, rays in muscle and marrow always describe longitudinal (pressure) waves. Conversely, cortical bone is considered as a solid. Rays in solids may represent longitudinal or shear waves. In the latter case, besides a power, a shear ray also has a polarisation vector, describing the oscillation direction of the particles. At an interface between two material types the incoming ray terminates and new reflected and transmitted (refracted) rays are generated. The directions of these reflected and transmitted rays are determined by Snell's law, the power distribution comes from the acoustic version of the Fresnel relations. At the liquid–solid interface like muscle–bone mode-conversion also takes place. The incoming ray from the liquid gives rise to three new rays: a reflected ray in the liquid, a transmitted (refracted) longitudinal ray in the solid and a transmitted (refracted) shear ray in the solid.

Figure 1 shows the power reflection/transmission coefficients at the muscle–bone interface, as calculated from the Fresnel relations, see the supplementary information and references given therein for the computation of the Fresnel relations (stacks.iop.org/PMB/61/1810/mmedia). Mode conversion also occurs in the opposite direction, i.e. at a solid–liquid interface like bone–muscle or bone–marrow. Shear waves in the solid may give rise to refracted longitudinal waves in the liquid. This also depends on the polarisation direction of the shear waves, see the supplementary information (stacks.iop.org/PMB/61/1810/mmedia). An overview of the ray-tracing process is shown in figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To compute the generated heat, each ray is divided into small segments of length $\Delta s$. The attenuation of acoustic power in a segment between distances s and $s+\Delta s$ from the initial point of the ray is given by

$$\begin{eqnarray}&&I(s)-I(s+\Delta s)=I(s)\left(1-{{\text{e}}^{-2\alpha \Delta s}}\right).\end{eqnarray}$$

Our model assumes that this power loss has occurred in the middle point of this segment (i.e. the point at distance $s+\frac{1}{2}\Delta s$). It is then assigned to the nearest grid point in a rectangular grid. The rays are terminated if their energy is less than a fraction (default 0.01) of the initial energy of a ray starting at one of the elements. In this way a three-dimensional table ${{\hat{Q}}_{ijk}}$ is computed, such that ${{\hat{Q}}_{ijk}}$ is the total power loss (in W cm⁻³) in the volume element around the grid point $\left({{x}_{i}},{{y}_{j}},{{z}_{k}}\right)$.

In fact, the power loss does not only give rise to heat production, other effects like scattering also play a role. This means that the actual heat production density Q_ijk is given by ${{Q}_{ijk}}=b\,{{\hat{Q}}_{ijk}}$, where b is the absorption fraction. The value of absorption fraction b depends on the tissue type and may be as low a 0.3 (Goss et al 1979). The ray tracer and the computation of the heat production density were implemented in MATLAB (MathWorks, Natick, MA).

Once the heat production density is known, the temperature T is found by solving the heat equation

$$\begin{eqnarray}&&\frac{\partial T\left(\mathbf{r},t\right)}{\partial t}=D\,\Delta T\left(\mathbf{r},t\right)+\frac{1}{\rho {{c}_{p}}}Q\left(\mathbf{r},t\right),\end{eqnarray}$$

where $Q\left(\mathbf{r},t\right)$ is the heat production density at position $\mathbf{r}$ and time t, and D is the heat diffusion coefficient, given by $D=\frac{{{k}_{\text{T}}}}{\rho \;{{c}_{p}}},$ where ${{k}_{\text{T}}}$ is the thermal conductivity, ρ is the density and c_p is the specific heat (at constant pressure), all depending on the material type. The value of $Q\left(\mathbf{r},t\right)$ is obtained from the heat production table Q_ijk (see section 2.1.1) and the sonication time. The sonication time is modelled by a step function, i.e. the heat source is described by Q_ijk under sonication and 0 otherwise. For an in vivo situation, the heat transport by perfusion has to be added to the heat source. We use von Neumann boundary conditions, which means that there is no heat flux through the boundaries of the system. The diffusion length in a time τ is approximately given by $\sqrt{4D\tau}$. For a simulation with a total duration of $\tau =100$ s and a diffusion coefficient in the muscle tissue of D  =  0.15 mm² s⁻¹, which gives a length of only 7.7 mm. Since the considered system is much larger than 7.7 mm, the temperature at the system boundaries will hardly change in 100 s, which justifies our assumption of no heat flux through the boundaries.

This heat equation can be solved by means of any standard finite-element software. We implemented our calculation in COMSOL Multiphysics (COMSOL, Inc., Burlington, MA).

The previously described approach, with a 1 mm radius target region for the rays, is used to model a 2 mm diameter sonication cell. In practice, larger regions are treated by superimposing 2 mm focal regions. A 4 mm diameter cell sonication is performed by superimposing 2 mm focal regions to cover the complete trajectory, as shown in figure 3. Since the (electronic) switching between the various 2 mm regions is very fast compared to the time scale of the temperature evolution, this process can be modelled by using one heat production density table Q, which is obtained by adding the tables for the individual 2 mm sub regions.

Zoom In Zoom Out Reset image size

This approach is not suitable for larger treatment cells. For instance, the 8 mm cell sonication in figure 3 is performed by dividing it into an inner and an outer circle of basic treatment cells (2 mm) according to the scheme presented in figure 3. This leads to two heat production densities: a table $Q_{ijk}^{(1)}$ corresponding to the inner trajectory 4 mm in diameter and a table $Q_{ijk}^{(2)}$ that corresponds to the outer trajectory 8 mm in diameter. When solving the bio-heat equation the actual heat source Q is, depending on the time, either Q⁽¹⁾ or Q⁽²⁾.

An alternative method to compute the heat production in one homogeneous material type is as follows. The pressure field of one transducer element is approximated by the far field approximation (Morse and Ingard 1986), which is a generally accepted approach. Then, by adding the pressure fields of all the transducer elements, the total pressure field and the corresponding heat production is obtained. Consider a circular transducer element with radius a, centered in the origin and directed towards the positive x axis, that emits acoustic waves with wave number $k=2\pi f/c$, where f is the frequency and c is the speed of sound. The far field approximation for the pressure field of this element at a position $\mathbf{r}=(x,y,z)$ is then given by

$$\begin{eqnarray}&&{{p}_{ffa}}(\mathbf{r})=\pi {{a}^{2}}{{S}_{0}}\frac{{{\text{e}}^{-\text{i}kr-\alpha r}}}{2\pi r}\,\frac{2{{J}_{1}}(ka\sin (\theta )}{ka\sin (\theta )}~, \end{eqnarray} \tag{ 1 }$$

where $r=\,|\mathbf{r}|$, θ is the angle between the vector $\mathbf{r}$ and the normal on the transducer element, J₁ is a Bessel function of the first kind, α is the (pressure) attenuation coefficient and S₀ is a measure for the strength of the element. The value of S₀ has been chosen such that the emitted power of each transducer element equals 1 watt, see section 2 in the supplementary information (stacks.iop.org/PMB/61/1810/mmedia).

Since the far field approximation assumes one homogeneous material type, this comparison is restricted to the case of a (theoretical) configuration with only one material type. For a theoretical comparison it is of no use to differentiate between absorption and attenuation. Hence, we assume here that for both the ray tracer and the far field approximation the absorption fraction b equals 1. In that case the heat flux and the heat production, as computed by the ray-tracing approach, are compared with the results from the far field approximation method. For the ray tracer the total heat flux through a coronal plane at distance x is given by

$$\begin{eqnarray}&&{{F}_{\text{R}}}(x)={{P}_{\text{init}}}-\int_{{{x}_{\text{min}}}}^{x}\left(\mathop{\int}^{}\mathop{\int}^{}\hat{Q}\left({{x}^{\prime}},y,z\right)\text{d}y\text{d}z\right)\text{d}{{x}^{\prime}}\end{eqnarray} \tag{ 2 }$$

where $\hat{Q}(x,y,z)$ is the heat production density computed by the ray tracer at position (x, y, z), ${{P}_{\text{init}}}$ is the initial acoustic power and the material with attenuation starts at position ${{x}_{\text{min}}}$. The integral is computed by summation over the corresponding part of the heat production table ${{\hat{Q}}_{ijk}}$, see section 2.1.1. Grid sizes of 2 mm and 0.5 mm were used for the ray tracer and the far field approximation, respectively.

Even if the total heat fluxes through the coronal planes of the two methods agree well, the distribution of the heat flux (and hence heat production) in the coronal planes may not be similar. Therefore, we also compared the integrated heat production over circles with radius R (in W/cm) in the coronal planes, see also figure 4. We used circles with radius R  =  8 cm, 3 cm, 1 cm or 0.4 cm. For both the ray-tracer and the far field approximation method the integrated heat production over parts of the coronal planes can be computed by integrating (summing up) over the corresponding part of the heat production densities Q and ${{Q}_{\text{F}}}$ respectively.

Zoom In Zoom Out Reset image size

For both methods the heat production was computed for $-10~\text{cm}\leqslant x\leqslant 2~\text{cm}$ and $-8~\text{cm}\leqslant y,z\leqslant 8~\text{cm}$, where the transducer is in the negative x direction, the focal distance is 12 cm and the natural focus point is in the origin, see also figure 2 for the used coordinates. The used parameter values are given in table 1. More details of the flux computation based on the far field approximation are given in section 1 of the supplementary information (stacks.iop.org/PMB/61/1810/mmedia).

Table 1. Parameters used in the simulation comparison.

f [MHz]	a [mm]	S0 [kg (m2 s2)−1]	ρ [kg m−3]	c [m s−1]	α [Np m−1]
1.4	3.3	$1.5007\centerdot {{10}^{7}}$	1040	1537	5.76

To assess the contribution of the shear waves to the total heat production in bone, two types of numerical simulations were performed. In the first simulations we compared the standard heat production in a bone to a hypothetical situation where the bone behaves acoustically as a liquid. In that case shear waves in bone do not exist and a much larger fraction of the acoustic energy is reflected.

For a bone a mesh was generated (by COMSOL) with a minimum element size of 0.1 cm and a maximum element size of 1 cm. The smaller elements were used in the focal region, where most of the heat is produced. In the mesh, two nested coaxial cylinders with radii of 1.0 cm and 2.0 cm, respectively, and height of 15 cm represents the marrow and the bone. (Thus, the thickness of the bone containing the marrow is 1 cm.) The bone cylinder is encased in a box with muscle properties with dimensions $\text{L}\times \text{W}\times \text{H}$ of $15\times 15\times 15~\text{c}{{\text{m}}^{3}}$. Finally, this box is nested within another box with water properties and dimensions $20\times 20\times 15~\text{c}{{\text{m}}^{3}}$. The distance from the transducer to the focus point and the centre of the mesh are 15 cm and 16.5 cm, respectively. The whole object can be rotated around the z-axis. The tilt of the object is measured by the angle γ with the z-axis. The configuration is similar to the one in figure 2. To be able to clearly show the tilt, the encasing muscle and water meshes are not shown in figure 2.

The rotation angle γ, i.e. the inclination of the mesh with respect to the z-axis, was varied in the range $-{{55}^{{}^\circ}}$ to 55° divided into steps of one degree. For each angle γ the transducer power was 30 W. The parameters for this experiment are given in table 2. To avoid computing the almost vanishing heat production in water and oil (due to the very small attenuation coefficients α) we have set the absorption fraction b to 0 for both materials. The heat production was calculated using a grid with a spacing of 0.1 cm only in the points residing in the bone cylinder. For each angle the ratio was computed between the total heat productions in the bone with and without shear waves. In the former case, to eliminate the shear waves, the bone was treated as a liquid.

Table 2. Parameters used in the simulations on shear waves and in the gel–bone comparison.

Material	c [m s−1]	ρ [kg m3]	cp [J (kg K)−1]	${{k}_{\text{T}}}$ [W (m K)−1]	α [Np cm−1]	b [ ]
Gel	1492	1030	3720	0.537	0.07	0.35a
Bone (longit)b	3736	2025	1300	0.487	1.9	0.3c
Bone (shear)b	1995	2025	1300	0.487	2.8	0.3c
Water	1523	997	—	—	$6.08\centerdot {{10}^{-4}}$	0
Oil	1380	1030	—	—	0	0

—The specific heat and thermal conductivity for water and oil are not needed, as solving the heat equation in these materials is unnecessary. ^aAverage of different values in Goss et al (1979). ^bNell and Meyers (2010). ^cPinton et al (2012).

In the second numerical experiment we investigated the effect of omitting shear waves on the resulting temperatures. The bone was defined as a cylinder 10 cm in height, with a 0.5 cm external radius and a 0.2 cm internal radius embedded in a box with $2\times 2$ cm² sides and 10 cm in height. The surrounding box was defined as a tissue mimicking gel. The transducer to bone distance and the focal distance were set to 13.5 cm. The bone geometry was meshed with a minimum element size of 0.1 cm. The computation was run for a fixed power of 120 W and two angles: a normal incidence and 50° incidence (close to total reflection limit). The other parameters were as in table 2, except the speed of sound in the tissue mimicking gel, which was 1537 m s⁻¹ in this simulation. The heat production was computed for each angle including and omitting shear waves mode-conversion. The heat equation solution was solved for each case considering 60 s sonication time followed by 60 s cooling time.

As a proof of concept, a bone phantom was designed matching a simple geometrical model. A pig femur was made into a cylindrical shape (figure 5(a)) with a cortex thickness of 0.7 cm. For temperature measurements inside the bone, three holes were drilled from the inside up to 3 mm below the cortex surface, to position the MR compatible optical temperature sensors (Neoptix T1 probes, Quebec City, Canada). The accuracy of the sensors was tested under different conditions (warm and hot water, free and within bone), using a laboratory thermometer as a reference.

Zoom In Zoom Out Reset image size

The bone was positioned in a beaker and then a tissue mimicking agar-silica gel was added to obtain the final setup, as shown in figure 5(b). The gel was prepared using 20 g agar and 30 g silica per litre of water, resulting in density, speed of sound and attenuation, as shown in table 1 (Partanen et al 2009). The phantom was positioned on a clinical HIFU table (Philips Sonalleve, Vantaa, Finland) in direct contact with the membrane on the tabletop, with water and additional gel pads to ensure acoustic coupling, in the centre of a 3T MRI scanner (Philips Achieva, Best, the Netherlands).

An MR-HIFU ablation experiment was performed using a protocol that is comparable to currently clinically used protocols for pain palliation procedures. For this protocol an 8 mm treatment cell (see section 2.1.3) was positioned 1 cm behind the bone surface. The HIFU transducer produced 110 W for 27 s, i.e. approximately 3 kJ energy. The incident angle was defined to be zero relative to the normal of the bone surface, as shown in figure 6. The focal region (figure 4(a)) was positioned along the central temperature sensor direction, thus measuring the temperature increase in the hottest region. The other sensors were positioned only 0.5 cm apart from each other to ensure that the heat diffusion could be monitored in time. The data from the sensors were recorded by a LabVIEW program (LabVIEW, National Instruments). MR-thermometry was performed using the PRFS technique (Denis De Senneville et al 2005, McDannold 2005, Rieke and Butts Pauly 2008, Köhler et al 2009), which is the state of the art for MR thermometry. Temperature maps were measured in 6 slices (see figure 6), and every 6 s a full temperature map for these slices was acquired (for more detail see table 3). MR-thermometry and the data from temperature sensors were afterwards compared with the simulations. The parameters used in the simulations are listed in table 2. To demonstrate the robustness of this model, three comparable experiments were performed using a 4 mm treatment cell, an acoustic power of 60 W, a treatment duration of 20 s and a frequency of 1.2 MHz. The exact experimental setup is described in the supplementary information under section 4 (stacks.iop.org/PMB/61/1810/mmedia).

Zoom In Zoom Out Reset image size

Table 3. MR thermometry—sequence parameters.

MR sequence M2D-EPI	TR/TE  =  38 / 20 ms
Spatial resolution	$1.42\times 1.42\times 4.08$ mm3
Temporal resolution	6 slices in 6.1 s