Ultrasound propagation modeling

Ultrasound are sound waves with frequencies higher than upper audible limit of human hearing, i.e. higher than 20 kHz. In the frequency range between 0.8 and 5 Mhz, the sound wave length is between 2 and 0.3 mm. Ultrasound is mainly used for diagnostic purposes as unlike light it can easily penetrate the tissue and is at those intensities completely harmless and is therefore used for diagnostics in pregnancies. Short wavelengths, enable us to create small regions of high intensity focused ultrasound at the focal point, which can be used for therapeutic purposes escoffre2016therapeutic such as for the tissue ablation. High intensity focused ultrasound has been already used in thermal ablation therapies escoffre2016therapeutic. In vital organs like brain it is particularly important to focus high intensity ultrasound beam at the small targeted structure in order to avoid unwanted tissue damage.

Based on our improvement of the of ultrasound responsiveness of the cells it would be feasible to introduce the genetic device into the target tissue and then selectively activate only cells by using focused ultrasound, where the sound pressure is sufficiently high. This would enable for example stimulation of neurons in the selected brain section to achieve noninvasive deep brain stimulation or stimulate cells that produce hormones.

For this purpose we designed a model of the ultrasound propagation to investigate the setup of the ultrasound device that we constructed, potentially using multiple ultrasound probes.

  Sound modeling with wave equation

The code for the ultrasound propagation model is written in C++ programming language with OpenGL library. The model was parallelized using OpenMP to take advantage of multi-core processor architecture and to achieve better performance. Source code is available on GitHub.

Sound is a vibration that propagates as a mechanical wave of pressure and displacement through the medium. The speed of its propagation depends on the properties of the medium. This means that speed of sound varies greatly between different types of tissue azhari2010basics. The speeds of sound in different tissue and in the air are listed in (2).

For a realistic simulation of the brain an MR image (3) was used from mollerpocket and the tissue segmentation was made by medical experts, our team members, students of medicine (Kosta, Samo and Nina) (4).

The probes in our model are capable of emitting sound waves at frequencies ranging from 0.5 Mhz to 5 Mhz at acoustic pressure from 10 kPa to 100 kPa. The focused probes that we used have focal length of 7 cm. Besides focused probes we also modeled unfocused probes. Both types have an aperture of 4.4 cm. In general propagation of sound from the focused probe can be approximated with an analytic function of a Gaussian beam (6), but this representation fails when we want to model wave interferences from multiple probes. For such cases wave equation model needs to be established.

Acoustic wave equation is a second-order linear partial differential equation that governs the propagation of acoustic waves through the medium. The general form of the equation is: \begin{equation} \frac{\partial^2 p}{\partial t^2} = c^2(\nabla^2 p) \end{equation} where $c$ is speed of sound in the medium, $p$ is acoustic pressure (the deviation from ambient pressure) and $\nabla^2$ is Laplacian operator (sum of second order partial derivatives with respect to each independent variable). In two dimensional Cartesian coordinate system we get the following equation: \begin{equation} \frac{\partial^2 p}{\partial t^2} = c^2(\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2}) \end{equation} While wave equation cannot be solved analytically, one can approximate it with finite difference method. Let's presume we want to approximate the solution on a square area. First we discretize the area to obtain the grid of evenly distributed points $(x_{i}, y_{j})$, $i,j = 1,...,n$. The Laplace operator can now be approximated: \begin{equation} \nabla^2 p_{i, j} \approx \frac{1}{h^2}(p_{i + 1, j} + p_{i-1, j} + p_{i, j + 1} + p_{i, j - 1} - 4p_{i,j}) \end{equation} where $h$ is Euclidean distance between two adjacent points. Since the frequencies between 0.8 and 5 Mhz correspond to the wavelength between 2 and 0.3 mm, $h$ does not need to be smaller than 0.3 mm. We now obtain a $n*n$ system of the second order ordinary differential equations: \begin{equation} \ddot{p}_{i, j} = \frac{c^2}{h^2}(p_{i + 1, j} + p_{i-1, j} + p_{i, j + 1} + p_{i, j - 1} - 4p_{i,j}) \end{equation} These equation can be simplified to first order differential equations with the introduction of new variable that represents the speed: $(v_{i, j} = \dot{p}_{i, j})$. The new system now has the following form: \begin{equation} \begin{aligned} \dot{p}_{i, j} = v_{i, j} \\ \dot{v}_{i, j} = \frac{c^2}{h^2}(p_{i + 1, j} + p_{i-1, j} + p_{i, j + 1} + p_{i, j - 1} - 4p_{i,j}) \end{aligned} \end{equation} The dispersion of waves is already included in wave equation, but absorption of sound waves in tissue is not included. However, absorption can be approximated with the introduction of the attenuation coefficient $k$. We fitted attenuation coefficient to experimental data sprawls1989ultrasound, $k = 16000$ for soft tissue and $k = 160000$ for bone tissue. The final equation now has a form: \begin{equation} \dot{v}_{i, j} = \frac{c^2}{h^2}(p_{i + 1, j} + p_{i-1, j} + p_{i, j + 1} + p_{i, j - 1} - 4p_{i,j}) - k*v_{i,j} \end{equation}

Finally we can solve given system with iterative numeric methods like Euler's method or Runge-Kutta method. The most widely known method of the Runge-Kutta family is RK4, which we also used in our model. While RK4 introduces some computational overhead it also significantly increases the numerical stability of the system in comparison to the Euler's method. The most basic steps used for the calculation of ultrasound propagation are displayed in (5).

We model a propagation of the ultrasound waves on a 500 x 500 grid, where distances between points are exactly 0.3 cm. This translates to a stimulated area of 15 by 15 cm. The model computes propagation of waves in 2 dimensions. While the computational complexity of calculating the solutions of wave equation in 3 dimensions is asymptotically much higher than in 2 dimensional space, it is unnecessary to demonstrate the ability to focus the ultrasound by modulation of the intensity, frequency and geometry of the ultrasound probes. However, in 3 dimensional space the intensity in the focal point should increase, since addition of sound waves is coming from all directions not just from left and right. Three or more probes could be combined to define the focal point in 3D. In our model it is also possible to observe the main sound characteristics such as the attenuation, reflection and occlusion. Probes are modeled as the set of finite elements. The pressure of every finite element on the probe is set according to sine function with the selected amplitude and frequency, thus every finite element acts like independent source of the acoustic pressure.

  Results

Simulation computed with our model (7) clearly shows that it is possible to target small regions of brain tissue. We can also observe the one dimensional cut of acoustic pressure, where focusing is even more evident (8). The probe on the right is slightly displaced due to asymmetry in tissue. This clearly shows that models need to be applied to calculate the optimal position of ultrasound probes. Such software will include automatic tissue segmentation, modeling of propagation of ultrasound waves and will assist experts in probe positioning. Our model is a step in that direction, since it can be used as a tool for positioning the probes. In the model that we generated focused and unfocused probes can be added or removed and positioned. Additionally probe's acoustic pressure and frequency can be set.

  Conclusion and outlook

We generated a model for the propagation of ultrasound waves through tissue. In our case a transversal cut of the middle head area was taken. Two probes were adjusted in such a manner that focusing on small area was obtained. With this model we can confirm the potential of the ultrasound as a platform for therapeutic applications, but for more advanced use the software tools needs to be further developed. That tool would include automatic tissue segmentation from patient medical images, computing propagation of ultrasound waves in 3D aimed to assist in experts in the probe positioning. The low intensity ultrasound image could be used to receive the feedback of the sonogram of the tissue for each probe that could be used to drive ultrasound stimulation for cell activation in the desired area. It is likely that the specific shape of the stimulated tissue could be selected rather than just a single focal point.

  References