Introduction

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. Unlike light ultrasound can easily penetrate the tissue and is at those intensities completely harmless. It is therefore used also in prenatal diagnostics where non-invasivness is of vital importance. Short wavelengths additionally enable us to create small regions of high intensity focused ultrasound waves at the focal point, which can be used for therapeutic purposes escoffre2016therapeutic for example 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 beams at the small targeted structure in order to avoid unwanted tissue damage.

Based on our improvement of the ultrasound responsiveness of the cells it would be feasible to introduce a genetic device into the target tissue and then selectively activate only specific cells by using focused ultrasound, causing a sufficient intensity of the acoustic pressure. This would for example enable the stimulation of neurons in the selected brain section to achieve non-invasive deep brain stimulation or to stimulate cells to produce different hormones.

For this purpose we designed a model of the ultrasound propagation. The model can guide us in the process of the ultrasound device calibration with the potential to use 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 also depends on the properties of the medium. This means that the speed of sound varies greatly between different types of tissues azhari2010basics. The speed of sound propagation in different tissue in comparison to its propagation through air is represented in (2).

In order to achieve a realistic simulation of the sound propagation through brain an MR image represented in 3 was used mollerpocket. Furthermore, we performed the tissue segmentation with the aid of medical experts, i.e. 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 of probes 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 squared area. First we need to discretize the area to obtain the grid of evenly distributed points $(x_{i}, y_{j})$, $i,j = 1,...,n$. The Laplacian operator can now be approximated as: \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 the first order differential equations with the introduction of a 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 the waves is already included in the wave equation. However, the absorption of sound waves in tissue is not included. The absorption can be approximated with the introduction of the attenuation coefficient $k$. We fitted attenuation coefficient to experimental data sprawls1989ultrasound to obtain, $k = 16000$ for soft tissue and $k = 160000$ for bone tissue. The final equation now has the 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 the given system with iterative numeric methods such as Euler's method or Runge-Kutta method. The most widely known method of the Runge-Kutta family is RK4, which was 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.