Modeling the interference between shear and longitudinal waves under high intensity focused ultrasound propagation in bone

High intensity focused ultrasound (HIFU) is a thermal technique emerged over the last years (ter Haar 2001). One important HIFU application is the treatment of various types of benign or malignant tumors (ter Haar 2001, Kim 2015). During HIFU cancers' therapy, the ultrasound (US) beam is focused on target tissues, such as the tumors, in the body (focal region), leading to rapid heating and subsequent coagulative necrosis (ter Haar and Coussios 2007, ter Haar 2011). HIFU treatments are usually performed under image-guidance using magnetic resonance imaging (MRI) or Ultrasound for planning and guidance. MRI-guidance has the advantage of measuring nearly real-time temperatures using MR-thermometry techniques such as the proton resonance frequency shift (PRFS) (McDannold 2005, Senneville et al 2005, Rieke and Pauly 2008, Köhler et al 2009). This method is suitable for water-rich tissues with a long transversal relaxation time (T2), and therefore works well for monitoring HIFU ablations of soft tissue tumors, such as uterine fibroids (Stewart et al 2003) and other emerging applications such as liver cancer (Al-Bataineh et al 2012, Leslie et al 2012, Aubry et al 2013). HIFU treatment of uterine fibroids focuses on alleviating the pain by reducing the tumor size (Vaezy et al 2000, Stewart et al 2006), in a safe and effective manner (Fischer et al 2015, Lee et al 2015), with acceptable adverse events (Wang et al 2014). MR-guided HIFU (MR-HIFU) has also been used as an alternative treatment strategy for palliative treatment of bone metastases (Li et al 2010, Hurwitz et al 2014, Huisman et al 2014, 2015, Chan et al 2016, Harding et al 2018). Bone metastases are associated with pain and reduced quality of life (Coleman 2006). General treatments focus on alleviating the pain, and aim at improving the patients' life conditions. MR-HIFU procedures have been shown to be technically feasible (Huisman et al 2014, Chan et al 2016), effective (Li et al 2010, Harding et al 2018) and safe with promising results (Huisman et al 2014). Clinical therapies with MR-HIFU have also been used recently for the treatment of osteoid osteoma (Napoli et al 2017, Sharma et al 2017, Yarmolenko et al 2018). Osteoid osteoma is a benign but painful neoplasm (Motamedi et al 2009) which occurs commonly in children and adolescents (Sharma et al 2017). MR-HIFU has been shown to be a valid alternative to the standard radiofrequency ablation techniques (Sharma et al 2017), while remaining feasible (Yarmolenko et al 2018), safe, and effective (Napoli et al 2017). Studies also state that HIFU treatments are a safe and feasible method for the treatment of malignant neoplasm osteosarcoma (Li et al 2009, 2010) and local unresectable recurrence of osteosarcoma (Yu et al 2015). However, the number of current treatments is small, and large-scale randomized control clinical studies are yet to be implemented (Yang et al 2018). MR-Thermometry (PRFS) is used during HIFU treatments to monitor whether a high enough temperature for ablation is achieved in the target region, while avoiding undesired heating of the surrounding healthy tissues. However, PRFS does not give the temperature information in bone due to low T2. Monitoring the temperature in bone is of utmost importance, given that the high ultrasound absorption in the cortex leads to an extremely rapid heating. Such rapid heating is undesired because it can cause unwanted side effects. Although innovative MR-based thermometry techniques for bone are already being investigated (Ramsay et al 2015a, 2015b), computer models can play a key role in the accurate prediction and monitoring of temperature. Modeling HIFU propagation in bone poses several challenges, due to the acoustic and thermal properties of the cortex (Goss et al 1979, Nell and Myers 2010). Moreover, at the interface of muscle and bone, all the reflected and transmitted ultrasound waves have to be considered. Particular attention must be paid to the shear wave mode conversion in bone. Since bone is a solid material (while muscle behaves as a liquid for ultrasound waves), both longitudinal and shear waves play a role in the heating process. Also, shear wave propagation cannot be neglected if the US waves hit the muscle-bone interface at an oblique incidence angle (Chan et al 1974, Fujii et al 1999, White et al 2006). Various studies including the shear waves HIFU propagation in bone have been performed (Fujii et al 1999, Lin et al 2000, Nell and Myers 2010, Treeby and Saratoon 2015, ten Eikelder et al 2016). Nevertheless, these methods tend to neglect the interference between longitudinal and shear waves or use simple transducer geometries. The latter enables the use of 2D-axisymmetric models. The use of 2D-axisymmetric models is sometimes necessary to avoid computer memory problems. However, it reduces the flexibility to adapt the model to transducers composed of multiple elements and to complex shaped propagation media (e.g. real bone CT data). In this paper, we present a computer ray tracing model for HIFU propagation, valid in soft tissues and bone. The model calculates the heat production density (in W cm⁻³) in the focal region and takes the interaction between longitudinal and shear waves into account. The model is able to deal with complex transducer geometries, and it can potentially handle propagation media with complex shapes, which makes it suitable for diverse medical applications. To investigate if our raytracer describes the HIFU propagation correctly, we compared the pressure in soft tissue and the displacement in bone computed by the raytracer with those obtained by solving the Helmholtz equation and the frequency-domain wave equation, in the case of one transducer element. Then, using a more realistic and complex configuration with multiple transducer elements (Philips Sonalleve system (Köhler et al 2009)) we investigated the role of the interference between longitudinal and shear waves, comparing the results given by two versions of the raytracer approach: one including the interference and the other neglecting it. In particular, we used two setups: the first one with only flat interfaces (lossless oil—soft tissue—bone) and the second one, more realistic, considering the bone as a cylinder with fat and marrow layers.

2.1. Modeling approach

2.2. Computation of the heat source: raytracer approach

In the raytracer approach, ultrasound waves are modeled as rays leaving the transducer element(s), with a certain frequency (harmonic oscillations are assumed). The orientation of the rays leaving each transducer element is given by a probability density calculated according to the far field approximation formula (Morse and Ingard 1968), which describes the pressure pattern in the far field zone generated by a circular transducer element (see figure 1). Each ray carries a power, a phase, a wave vector (parallel to the propagation direction of the rays) and, in case of shear rays in bone, also a polarization direction. The spatial distribution of the rays emitted from each transducer element is chosen such that the average intensity emitted by the rays corresponds to the far field distribution of the transducer element. The power of the ray $\newcommand{\e}{{\rm e}} W(\ell)$, after a distance $\newcommand{\e}{{\rm e}} \ell$ in a material with an attenuation coefficient for pressure $\alpha$, is computed as $\newcommand{\e}{{\rm e}} W(\ell)=W(0)e^{-2 \alpha \ell}$. At each interface, the transmission and reflection coefficients are computed, and all the reflected/refracted rays are generated, taking a possible phase shift into account. Rays showing the same pattern (i.e. which are generated by the same transducer element, and with almost the same physical path to an observation point $\boldsymbol{r}_0$) describe waves propagating from the transducer element to the observation point. Rays with the same pattern have the same reflection/refraction sequence. To find the contribution of the waves with this pattern to the total pressure or displacement, the intensity and phase from the rays in the pattern have been computed. That results in an expression for the pressure or displacement generated by the considered pattern. Since the considered wave equations are linear, the pressures (in soft tissue) and displacements (in solids) corresponding to the various patterns can be added to obtain the total pressure or displacement. Finally, the total pressure or total displacement is used to compute the total power loss due to attenuation, which results in heat production. In the following subsections more details of the various steps are given.

Zoom In Zoom Out Reset image size

2.2.1. Intensity for soft tissue and bone

Consider a point $\boldsymbol{r}_0$ in soft tissue or in bone, and a volume (cube) around it. When a ray of a given pattern $g$ intersects the cube, the part of the ray inside the cube is described as an infinitesimally small cylinder with a circular base area $\Delta A_l$ and a height $s_l$ (see figure 2). The intensity associated to the ray in the cylinder is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_l= \frac{W_l}{\Delta A_l}, \label{int} \nonumber \end{align} \tag{ 1 }$$

where $W_l$ is the power of the ray. Equation (1) holds only in the cylinder with circular base area $\Delta A_l$. Considering that the volume of the cylinder enveloping the ray is given by $\Delta A_l s_l$, the contribution of a ray to the total average intensity in the cube is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_l \frac{\Delta A_l s_l}{\Delta V} = \frac{W_l}{\Delta A_l} \frac{\Delta A_l s_l}{\Delta V}=\frac{s_l}{\Delta V} W_l, \label{int3} \nonumber \end{align} \tag{ 2 }$$

where $\Delta V$ is the volume of the cube. The total intensity in the cube around $\boldsymbol{r}_0$, for all the incoming rays with the same pattern, is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle I_{(g)}= \sum_{l=1}^{N_r} I_l , \label{int44} \nonumber \end{align} \tag{ 3 }$$

where $N_r$ is the total number of rays of pattern $g$ intersecting the small cube around $\boldsymbol{r}_0$.

Zoom In Zoom Out Reset image size

2.2.2. Pressure in soft tissue

For soft tissue, the (complex) pressure ${p}_{(g)}$ in a point $\boldsymbol{r}_0$ is related to intensity given in formula (3) as:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {p}_{(g)}= \sqrt{2ZI_{(g)}}e^{i \phi_{(g)}} , \nonumber \end{align*}$$

where $Z$ is the acoustic impedance and $\phi_{(g)}$ is the averaged phase in the small cube around $\boldsymbol{r}_0$ calculated as:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \phi_{(g)}= \frac{\displaystyle \sum \nolimits_{l=1}^{N_r} \phi_l}{N_r} , \nonumber \end{align*}$$

where $\phi_l$ is the phase associated to the longitudinal $l$th ray.

2.2.3. Displacement in bone

The intensity $I_L$ for longitudinal waves of the same pattern $g$ in bone is related to the absolute amplitude $A$ of the displacement associated to longitudinal waves by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {I_L}_{(g)}= \frac{1}{2}\omega~A^2( (\lambda + 2\mu)k_L + (\xi +4\eta/3)\omega \alpha_L),\label{uff1} \nonumber \end{align} \tag{ 4 }$$

where $\omega$ is the angular frequency. $k_L$, $\alpha_L$ and $c_L$ are the plane wave parameters related to the the Lamé coefficients $\lambda$ and $\mu$, to the first and second viscosities $\xi$ and $\newcommand{\e}{{\rm e}} \eta$ and to the material density $\rho$ (see section 1 in supplementary information (stacks.iop.org/PMB/63/235024/mmedia)).

In the same fashion, the intensity $I_S$ for shear waves of the same pattern $h$ in bone is related to the absolute amplitude $B$ of the displacement associated to shear waves by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {I_S}_{(h)}= \frac{1}{2}~B^2 (\mu \omega k_S + \eta \omega^2 \alpha_S),\label{uff2} \nonumber \end{align} \tag{ 5 }$$

where $k_S$, $\alpha_S$ and $c_S$ are the plane wave parameters related to $\lambda$, $\mu$, $\xi$, $\newcommand{\e}{{\rm e}} \eta$ and $\rho$ (see section 2 in the supplementary information). An explicit derivation of relations (4) and (5) can be found in sections 4 and 5 of the supplementary information. The absolute amplitudes $A$ and $B$ can be easily found from equations (4) and (5), respectively.

Then, using the averaged phases for longitudinal $\phi_L$ and shear waves $\phi_S$ that are also obtained from the rays in the considered pattern, the complex displacements around $\boldsymbol{r}_0$ can be written as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{\boldsymbol{u}}_{L(g)}&= A~e^{i \phi_L} \boldsymbol{\tilde{k}}_L, \nonumber \\ \hat{\boldsymbol{u}}_{S(h)}&= B~e^{i \phi_S} \boldsymbol{a}, \label{disp} \nonumber \end{align} \tag{ 6 }$$

where $\boldsymbol{\tilde{k}}_L$ is the averaged and normalized wave vector for longitudinal waves with the same pattern $g$ and $\boldsymbol{a}$ is the averaged and normalized polarization vector for shear rays with the same pattern $h$.

2.2.4. Calculation of the total heat production density

The total pressure $p$ (in soft tissue) and the total displacement $\hat{\boldsymbol{u}}$ (in bone) are obtained by adding the complex pressures and displacements of the various patterns. Adding complex values ensures that the interference is included.

Finally, the total power loss can be computed using:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Q_bn_0} Q=(\xi + \eta/3)(\sum_l \dot{u}_{l,l})^2+\eta \sum_{l,j}{(\dot{u}_{l,j})}^2+ \eta \sum_{l,j}(\dot{u}_{l,j}\dot{u}_{j,l}-\dot{u}_{l,l}\dot{u}_{j,j}), \nonumber \end{align} \tag{ 7 }$$

where $\dot{u}_{l,j}$ represents the differentiation of the $l$th component of the velocity with respect to coordinate $j$. Formula (7) is still time-dependent, and therefore its time average $\langle Q \rangle$ has been computed. More information about the derivation of formula (7) can be found in section 3 of the supplementary information. The spatial derivatives of the displacements in expression (7) are explained in sections 6 and 7 of the supplementary information. It is worth mentioning that the time average of equation (7) for soft tissue is reduced to the standard formula:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \left<Q \right>= \alpha \frac{{\left\Vert\,p \right\Vert }^2}{Z} , \nonumber \end{align*}$$

where $\alpha$ is the attenuation coefficient for soft tissue. The heat production density, which gives rise to heat production, is found as:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \hat{Q}= b_a \left<Q \right> , \nonumber \end{align*}$$

where $b_a$ is the absorption fraction, defined as the actual fraction of power loss converted into heat, which depends on the tissue type. In fact, not all the power loss is responsible for the temperature increase, but also other effects like scattering take place. Typical values of $b_a$ for soft tissue and bone can be found in table 1.

Table 1. Parameters used in simulations.

Material	$c$ (m s−1)	$\rho$ (kg m−3)	$\alpha$ (Np cm−1)	$b_a$ (—)	$c_p$ (J (Kg K))$^{-1}$	$k_T$ (W (m K)$^{-1}$)
Soft tissue—muscle—marrow	1500	1050	0.07	0.35a	3720	0.537
Fat	1450	909	0.034	0.35a	—	—
Bone longitb	3736	2025	1.9	0.3c	1300	0.487
Bone shearb	1995	2025	2.8	0.3c	1300	0.487
Lossless oil	1380	1030	0	—	—	—

^aAverage of different values in Goss et al (1979). ^bNell and Myers (2010). ^cPinton et al (2012).

2.3. Comparison with a finite element approach with only flat interfaces

To investigate if the raytracer approach describes the US propagation in the correct way, a comparison between the raytracer and the solutions of the Helmholtz equation in soft tissue and the frequency-domain wave equation in bone has been done. In the soft tissue portion, the Helmholtz equation with the total pressure $p$ as unknown is written as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \triangle p - k^{2} p =0, \label{Helmholtz2} \nonumber \end{align} \tag{ 8 }$$

where $k$ is the (complex) wave number for soft tissue defined as $k= \frac{\omega}{c} - i\alpha$, where $c$ and $\alpha$ are the speed of sound and attenuation for the soft tissue. In bone, the equations to solve are:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} - \omega^2 \rho \boldsymbol{u}=\nabla \cdotp \boldsymbol{s} \label{momentum} \nonumber \\ {s}_{lj} = \sum \limits_{m,n} {C}_{ljmn}~\epsilon_{mn} \nonumber \\ \boldsymbol{\epsilon}= \frac{1}{2} ( \nabla \boldsymbol{u} + (\nabla \boldsymbol{u})^T). \end{array} \right. \nonumber \end{align} \tag{ 9 }$$

where $\newcommand{\e}{{\rm e}} \boldsymbol{\epsilon}$ is the strain tensor, $\rho$ is the density, ${s}_{lj}$ are the components of the stress measure at each point, $\boldsymbol{u} (u,v,w)$ is the displacement vector, and ${C}_{ljmn}$ are the component of the fourth-order stiffness or elasticity tensor, which are dependent on the Young modulus $E$ and the second Lamé coefficient $\mu$, also known as shear modulus. In bone, the attenuation is taken into account by including damping:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle E = E' + i~\eta_e E',\nonumber \\ \mu = \mu' + i~\eta_g \mu', \nonumber \end{align*}$$

where $\newcommand{\e}{{\rm e}} \eta_e$ is the loss factor for Young modulus $E$, $E'$ is the real part of $E$, $\newcommand{\e}{{\rm e}} \eta_g$ is the loss factor for shear modulus $\mu$ and $\mu'$ is the real part of $\mu$.

To define the same attenuation for the raytracer model and for the solution of equation (9), a relation between the attenuation coeffients $\alpha_L$ and $\alpha_S$ and the two loss factors $\newcommand{\e}{{\rm e}} \eta_e$ and $\newcommand{\e}{{\rm e}} \eta_g$ has been derived (see section 8 in the supplementary information).

Equation (9) have been implemented for one circular transducer element with radius $2$ mm, using a 2D-axisymmetric model in COMSOL Multiphysics (Version 5.0, Burlington, MA, 2014). The transducer element is immersed in soft tissue, and its center is in (0,0,0) cm. The ultrasound propagation is along the $x$-axis and the flat interface with bone is located at $x = 1.6$ cm. A schematic representation of the simulation setup and of the boundary conditions used is given in figure 3. The comparison with the raytracer has been done in terms of the pressure in the soft tissue and in terms of the displacement field in the bone. The frequency used is 1.2 MHz and the initial power is 1 W. The soft tissue and bone parameters are shown in table 1. A comparison between the raytracer and the solution of the frequency-domain wave equation has also been implemented by modeling the bone as a cylinder (see section 9 in the supplementary information).

Zoom In Zoom Out Reset image size

2.4. Computation of the temperature evolution

When the total heat production density is known, the heat equation can be solved to find the temperature evolution in the desired volume. The heat equation is formulated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\tr}{{\rm tr}} \displaystyle \dot{T}(\boldsymbol{r},t)= D \triangle T(\boldsymbol{r}, t) + \frac{1}{\rho c_p} \hat{Q}(\boldsymbol{r},t), \label{temp} \nonumber \end{align} \tag{ 10 }$$

where $T$ is the temperature at position $\boldsymbol{r}$ and time $t$, $D$ is the diffusion coefficient defined as $D= \frac{K_T}{\rho c_p}$ with $K_T$ the thermal conductivity, $\rho$ the density and $c_p$ the specific heat at a constant pressure. The value of $\hat{Q}(\boldsymbol{r},t)$ has been modeled by a step function, i.e. it is equal to the heat production density computed by the raytracer in the point $\boldsymbol{r}$ under sonication, and zero otherwise. We assumed that there is no heat flux through the boundaries of the system (Von Neumann boundary conditions) both in soft tissue and in bone. To obtain the temperature rise during HIFU, a finite-element software package has been used (COMSOL Multiphysics, Version 5.0, Burlington, MA, 2014).

2.5. Contribution of shear and longitudinal waves interference

To investigate the role of the interference between longitudinal and shear waves during a HIFU procedure, the geometry of the phased array transducer in the Philips Sonalleve system has been reproduced by the raytracer model. The system is composed of 256 transducer elements, positioned in a spherical shell and immersed in a lossless material. Each circular element has a radius of 3.5 mm. The system operates at 1.2 MHz, and it creates a cigar-shaped focal region with a dimension of approximately 2 mm width and 7 mm length. Two numerical simulations for two different setups have been performed (see figure 4). The first setup considers flat interfaces only, while the second one models a more realistic case, including a fat and a muscle layer, as well as a cylindrical bone. The used working frequency is 1.2 MHz, the initial total power is set at 25 W for the first setup and 40 W for the second one. The material parameters are given in table 1. In the first simulation, the total heat production density is computed in bone as a sum of $Q_L$ and $Q_S$, where $Q_L$ is the total power loss when only longitudinal waves are present and $Q_S$ is the total power loss calculated when only shear waves are present. In this case, the interference between longitudinal-longitudinal and shear-shear and the contribution of shear waves are still taken into account, while the interference between shear and longitudinal waves is neglected. In the second simulation, the total power loss in bone is computed using the time averaged version of formula (7), containing the interference information. A comparison of the total power loss integrated over different volumes (in W) has been done using the first setup for the two simulations. Each volume is represented by a 3D half ellipsoid with semi-axes in $y$- and $z$-direction ($b = c = 0.4$ cm), and a length of the remaining semi-axis ($a$ in the $x$-direction) between 0 cm and 0.7 cm (see figure 5). Further, the temperature rise has been computed for both setups, with the heat source applied for 10 s. Both the total power loss and the temperature rise are computed on a grid with a maximum spacing of 0.29 mm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Comparison with a finite element approach

Figure 6 shows the pressure trends computed by solving the Helmholtz equation for soft tissue and the raytracer approach for the configuration in figure 3. The pressure was computed along the central line from the centre of the transducer element in the $x$-direction. As the raytracer model was based on the far field approximation approach, it is considered to give a correct description of the ultrasound propagation only in the far field zone. The near field distance in the US propagation direction, for a circular transducer element of diameter $d$, was calculated as $d^2 / (4 w)$, where $w$ is the wavelength. The near field distance with $d= 4$ mm and $w= 1.25$ mm was found to be at 3.2 mm. For $x$-values smaller than 3.2 mm, the behavior of the pressure is not shown. In the far field zone, the two models described an output with a virtually identical oscillation and similar pressure values.

Zoom In Zoom Out Reset image size

In figure 7, the $x$- and $y$-components ($d_1$ and $d_2$) of the absolute value of the displacement in the $xy$-plane of bone given by the raytracer and by solving equation (9) for the configuration in figure 3 are shown. The displacement in the $z$-direction has been omitted because it was zero everywhere for both the methods. The peak values for $d_1$ were 8.74 · 10⁻⁷ cm and 9.02 · 10⁻⁷ cm, and for $d_2$ were 2.98 · 10⁻⁷ cm and 3.02 · 10⁻⁷ cm for equation (9) and the raytracer respectively. The two models predicted comparable displacements with regard to the peak values and to the field distribution in the planes, e.g. the low values assumed by $d_2$ in the central region of the $xy$-plane were visible in both the methods.

Zoom In Zoom Out Reset image size

3.2. Contribution of shear and longitudinal waves interference

In figure 8, $Q_{v}$, the total predicted power loss integrated over different volumes for the raytracer with ($R_{int}$) and without ($R_{sum}$) interference between longitudinal and shear waves, is shown. Computing $Q_v$ with $R_{sum}$ led to an over-estimation of the integrated power loss of about 75%.

Zoom In Zoom Out Reset image size

In figures 9 and 11, the predicted temperatures in the $xy$-plane at time 11 s for $R_{int}$ and $R_{sum}$ for the two setups are shown. When using the first setup, the maximum temperature increase in bone predicted by $R_{sum}$ was approximately 3.4 times higher than the one given by $R_{int}$. The over-estimation of $R_{sum}$ was also visible when using the second setup, where the temperature predicted by $R_{sum}$ was around 4.1 times higher than the one given by $R_{int}$. Higher values of temperature for $R_{sum}$ were also observed when investigating the temperature trend over time at point ($-1$,0,0) cm (see figures 10 and 12) for the two setups.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We presented a raytracer model that describes the HIFU propagation in soft tissue and bone. Our aim was to investigate the interaction between longitudinal and shear waves in bone. The method was compared with the solution of the Helmholtz equation and the frequency-domain wave equation, in the case of one transducer element. Figure 6 shows that the pressure field predicted by the raytracer and by solving the Helmholtz equation agree well in the far field zone. At a distance within the near field length, the raytracer method is not applicable, because the rays are generated according to the far field approximation pattern. As our model assumes systems that focus on tumors in the far field area, e.g. the Philips Sonnalleve, it is logical to focus on the correct description of the US propagation only in the far field zone. Figure 7 shows a pattern in the displacement field inside the bone, both for the raytracer and the wave equation (e.g. higher displacement in the $x$- direction close to the interface in the central portion, and in the $y$- direction close to the interface in the upper and lower portions). The results compare very well with only minor differences in the values of the displacement in the $x$-direction (see figures 7(a) and (b)). These minor differences could be due to the small approximation error inherent to the discrete nature of the raytracer approach (i.e. the finite number of rays used in the simulations).

The inclusion of shear propagation is crucial when the angle of incidence ($\theta_i$) of US waves on the interface is oblique (Chan et al 1974, Fujii et al 1999, White et al 2006). Moreover, when $\theta_i$ is between 21° and 52°, the energy transferred to the bone by shear waves is more than the energy transferred by the longitudinal waves in the same interval (ten Eikelder et al 2016). In the Philips Sonalleve system, the angle between the normal to the outermost transducer elements and the normal to the coronal plane is more than 32°. Therefore, in our system, the shear wave propagation cannot be neglected. Furthermore, figures 8–12 show that it is important to include in the model both shear wave propagation and interference (constructive/destructive) between shear and longitudinal waves. Figure 8 in particular shows that not including the interference between the two waves leads to an over-estimation of the integrated power loss, causing an over-estimation of the temperature increase (figures 9–12). Not including the interference between shear and longitudinal waves gives an over-estimation of the temperature increase of about 3.4 times when using the first setup and of about 4.1 times when using the second setup. This slight difference can be due to the presence of the curved surface in the second setup, which increases the presence of shear waves and consequently the interference between shear and longitudinal waves grows in importance. The general over-estimation given by not considering the interaction between shear and longitudinal waves may be explained by the different phases of longitudinal and shear waves, which give rise to destructive interference. For US waves originating from the very same transducer element, the interference between longitudinal and shear waves is partially caused by the different sound velocities that these two kinds of waves have. At the interface of soft tissue and bone, the generated refracted shear waves show a phase shift with respect to the incoming wave, which is not observed for longitudinal waves. When the US waves originate from different transducer element, also the relative position of both transducer elements plays a role. For shear waves, the interference also occurs because of a difference in polarization directions. Hence, for HIFU propagation in bone, the interaction between shear and longitudinal waves plays an important role for the temperature evaluation in the focal region.

Several studies about the role of shear waves in HIFU propagation on the temperature output for bone have been performed (Fujii et al 1999, Lin et al 2000, Nell and Myers 2010, Treeby and Saratoon 2015, ten Eikelder et al 2016). Fujii et al (1999) did not include longitudinal and shear waves interaction in bone, as they focused on the temperature evolution in the portion of muscle close to the interface with bone. Lin et al (2000), Treeby and Saratoon (2015) and ten Eikelder et al (2016) calculated the total heat production density in bone by summing the contribution of longitudinal waves and shear waves without including interference between the two waves. Nell and Myers (2010) considered bone as a solid-elastic material and they used a finite-element commercial software to calculate the temperature, including the interference between longitudinal and shear waves. However, the used transducer geometry is far less complex than more realistic transducer geometries, like the Philips Sonalleve transducer geometry that we used in this study. In fact, it is composed of one spherically focused transducer with a diameter of 8 cm, while the Philips Sonalleve system comprises 256 transducer elements (with 12 cm radius of curvature, 13 cm aperture). Nell and Myers (2010) used a 2D-axisymmetric model, not suitable for realistic transducer geometries, which aim to have as little geometric symmetry as possible to avoid overheating in the near-field (Ramaekers et al 2017) as well as to reduce the presence of grating lobes (Pompei and Wooh 2002). This type of models is also not adequate for complex shaped propagation media (such as real bone). In summary, the raytracer approach is flexible and it can easily deal with various transducer configurations. Moreover, it can handle 3D propagation geometries composed of multiples material layers, without becoming dramatically more memory intensive. This makes the use of the raytracer preferable to software packages for US propagation (e.g. COMSOL), as they often show a high memory consumption when treating 3D geometries. Our method has also the potential of handling more complex shapes (e.g. CT data, complex geometry described by meshes), which are able to describe real bones. In the future, the experimental validation of the raytracer approach using CT data should be implemented. Currently, considering the configuration with 256 transducer elements, the MATLAB version of the code takes around 50 min to generate the power production (in 56 000 points) when a cylindrical bone is present, and around 20 min when only soft tissues are considered. Therefore, we are currently working on the optimisation of the code, to reduce the computation time. This may lead to the application of the model for real bone lesions of patients, extending its usage to a clinical environment (i.e. treatment planning).