2.1

HIFU-exposed tissue phantom model and governing equations

A schematic drawing of the model used for simulating the acoustic pressure in the tissue exposed to HIFU is shown in Fig. 1. 3D geometry of the model was built on a Cartesian space coordinate. A tissue-mimicking phantom was constructed, consisting of skin, fat, muscle, and blood layers, with a focus on the peripheral artery disease such as iliofemoral artery occlusion. Referring to the clinical usage of the HIFU with insonation frequencies of 0.5-2.0 MHz for the treatment of thrombo-occlusive disease, single element HIFU transducers of H-104, H-101, H-199, and H-106 bowl types (Sonic Concepts, Bothwell, WA, USA) were utilized as a HIFU source. These transducers have central frequencies of 0.5, 1.1, 1.5 and 2.0 MHz, respectively. An acoustic wave generated by a HIFU source, with an aperture radius of 32 mm and curvature radius of 63.2 mm, propagates through water and tissue, and is focused on the blood vessel, which is prolonged alongside x-axis at the focus (z = f_c = 63.2 mm). The thicknesses of the skin, fat, and muscle layers were set to be 2, 3, and 31.8 mm, respectively (the surface of the tissue was located in z = 53.2 mm), and the diameter of the blood vessel, d_v, was 2 mm. These dimensions were fixed in the entire simulation unless otherwise to account for individual size effects. Considering the plane symmetry of the domain on the xoz and yoz plane, we have only utilized one quarter of the total size in simulation to reduce computational costs.

Figure 1.

Schematic diagram of HIFU-exposed tissue phantom for simulation.

The nonlinear propagation of HIFU is described by full-wave Westervelt equation, which considered the influence of diffraction, absorption, and nonlinearity, given by, [4]

where p denotes the acoustic pressure, t is time, ∇² = ∂²/∂x² + ∂²/∂y² + ∂²/∂z² is the Laplacian operator, c and ρ are the speed of sound and mass density of the medium, β and δ represent the nonlinear coefficient and sound diffusivity of the medium, respectively. The sound diffusivity is calculated by,

where α means the absorption coefficient of the sound and ω is an angular frequency of the acoustic wave. The absorption coefficient depends on the acoustic frequency, f, as follows,

where α0 is the absorption coefficient at f₀ = 1 MHz and the exponent γ in the water and in the tissue is 2 and 1.1, respectively [7]. To prevent spurious reflections from the boundaries of the computational domain, absorbing boundary conditions were accounted by assuming a perfectly matched layer (PML) outside of the model (see Fig. 1).

2.2

Simulation and its convergence test

The nonlinear Eq. (1) was solved for the acoustic media using the finite-difference time-domain (FDTD) method [4]. To model the pressure on a surface of focused circular transducer, planar pressure source approximation [4] has been applied at z = 0. Zero initial condition of the pressure and its derivative with respect to time (i.e., p(x, y, z)=0, ∂p / ∂t = 0) has been set throughout whole domain of the physical body. The absorbing boundary condition in the PML has been implemented following the steps described by Gupta & Srivastava [4]. While, symmetrical conditions, i.e., p(x, -y, z, t) = p(x, y, z, t), p(-x, y, z, t) = p(x, y, z, t), have been set on the xoz and yoz plane, respectively. The pressure amplitude of the HIFU source, p₀, was obtained by Eq. (4),

where A represents the area of the transducer surface, P_act = ηP_elect is the acoustic power of the transducer, Pelect is the electrical power supplied to the transducer. The nominal conversion efficiency of the electrical to acoustic power, η, is 0.85 in H-series transducers, measured by manufacturers. Acoustic properties of the medium used in this simulation are listed in Table 1.

Table 1.

The acoustic properties of the medium [4]

To test convergence of the solution depending on the numerical parameters (discretization steps of space and time, Δx, Δy, Δz and Δt), simulation for the water medium has primarily been performed varying with the parameter sizes in strict. Figure 2 shows the pressure variations along the z-axis and x-axis for different discretization step sizes in space, and according to time for different time step size. When Δt was set to T/100 (T is period of wave), the pressure distribution exhibited a satisfactory agreement with variation of Δz from λ/10 to λ/12 [λ is wavelength, see Fig. 2(a)], and of Δx and Δy from λ/5 to λ/7 [see Fig. 2(b)]. The pressure alternation in the cross-section to wave propagation is not such serious compared to that of axial direction (z-axis). Therefore, even if the grid resolution in x-axis and y-axis was not higher than in z-axis, good convergence of the solution could be accomplished. In addition, when Δt changed from T/50 to T/120 (Δx=Δy=λ/5, Δz=λ/10), almost no differences were found in magnitude and shape of curves [Fig. 2(c)]. Through such convergence test of the solution, optimized discretization step sizes, Δx, Δy, Δz, and Δt were set to λ/5, λ/5, λ/10, and T/100 respectively. The simulations were performed in MATLAB (R2018a) on a desktop PC with Intel (R) Core (TM) i7-7700 CPU of 3.6 GHz and 64 GB RAM. Depending on the frequency, it required approximately 15 minutes to 4 days to complete one simulation.

Figure 2.

Convergence test of the solution conducted at f = 2 MHz and pelect = 150 W, comparing with instant pressure distribution (a) along the z-axis for different Δz (Δx=Δy=λ/5, Δt=T/100), and (b) along the x-axis for different Δx and Δy (Δz=λ/10, Δt=T/100), and also (c) with pressure waveform for different Δt (Δx=Δy=λ/5, Δz=λ/10).

2.3

Simulation under different conditions

After completing the convergence test, the simulation was performed under the condition of P_elect = 60 W for the tissue model depicted in Fig. 1, to visualize the distribution of the pressure field in the tissue. Subsequently, the simulation was repeated with P_elect varying from 30 W to 150 W, to evaluate the peak pressure in the blood region based on the insonation power and frequency.

To account for different lesions and individuals, the simulation was also conducted under various tissue conditions by varying the depth of blood vessels, and thickness of the tissue (skin and fat) layers, respectively, while keeping the power fixed at P_elect = 120 W. The depth of the blood vessel (from the surface of the tissue to a vessel center) ranged from 10 mm to 30 mm, with the remaining mediums set as the same values as in Fig. 1. Lastly, the thicknesses of the skin and fat layers varied from 1-5 mm and 2-7 mm, respectively, toward the transducer surface, while keeping the muscle thickness at 31.8 mm and the vessel location at z = f_c.

3.1

Pressure distribution in the tissue

Figure 3 shows the pressure field distribution in the tissue obtained at P_elect = 60 W. In the focal region, the peak positive pressure, p⁺, is greater than the peak negative pressure, p^-, and this difference becomes more pronounced with increasing frequency. Interestingly, a difference in the acoustic properties between media caused multiple reflections with the layered media, showing the tooth shape of the curve [see Fig. 3(a)]. It is expected that serious reflections are generated in the interfaces of the skin-fat and fat-muscle layers, since they have a relatively large reflectance (|V_p| = 0.14, 0.13 when the degree of the incident angle of the beam, θ, is 30°). Similarity of the physical properties (density and sound speed) of the muscle and blood lead to a weak reflection (|V_p|=0.0038 at θ=30°), thereby not observing noticeable fluctuation in the interfaces between media. Hence, in the cross-sections [xoy plane in Fig. 3(b)] involving the muscle and blood region, the pressure distribution exhibits a smooth appearance, without showing any leap patterns. Although ultrasound propagates through multiple media layers, the shape and location of the focus showed no serious changes from that appeared in the water at all frequencies, i. e., the geometric focus on the blood vessel did not shifted, unlike demonstrated in the case of passing through cerebral tissue covered with skull [8]. Increasing the frequency improved the focusing characteristics with shorter wavelength, thus leading to an increment of the peak pressure.

Figure 3.

(a) Peak pressure contour and curve on yoz plane and along the z-axis, respectively, (b) peak pressure contour and curve on xoy plane and along the x-axis at z = fc, respectively, obtained at Pelect = 60 W. The vertical and horizontal dotted lines represent the interfaces between adjacent mediums.

3.2

Pressure changes with different HIFU powers and frequencies

Figure 4 shows the results of the peak pressure in the blood region (z = fc) simulated at different HIFU powers and frequencies. Both the peak positive pressure, p⁺, and peak negative pressure, p^-, increase as the power of HIFU transducer increases. The increment of p⁺ is more dramatic compared to that of p^-, implying that the nonlinearity of HIFU is strengthened at a higher intensity. Moreover, as demonstrated in Fig. 3, the magnitude and increment of the peak pressure also increase as the frequency increases.

Figure 4.

Peak pressure in the blood region (z = fc) versus an electrical power, P_elect, and frequency, f, supplied to HIFU transducer

3.3

Effect of depth of blood vessel

Figure 5 is the simulation results for the peak pressure in the blood region, where the depth was changed from 10 mm to 30 mm. The simulation was performed at three different HIFU outputs, i.e., P_elect = 60, 90, and 120 W. It is observed that the deeper the depth of the focus, the smaller the pressure becomes. The results demonstrate that the attenuation of the pressure varies in relation to both the insonation power and the frequency. In the case of p⁺ value, the attenuation is relatively small when f ≤ 1.1 MHz, whereas for f > 1.1 MHz, the attenuation strengthens with an increase of the HIFU output. On the other hand, for p^- value, the attenuation is not so remarkable for all frequencies and HIFU outputs. For instance, at Pelelct = 120 W and f = 2 MHz, the acoustic attenuation for p⁺ is 29.0 Np/m, while for the p^-, it is approximately half of that with p⁺, measuring 15.2 Np/m when calculated by -1/Lln(p/p0) in Fig. 5(c), where p0 and p are the acoustic pressure at z = z0 and z = z0 + L.

Figure 5.

Peak pressure changes with the depth of the blood vessel at (a) P_elect = 60 W, (b) 90 W, and (c) 120 W.

3.4

Effect of tissue thickness

Figure 6 represents the result of the simulation for the peak pressure inside the blood vessel (z = f_c) as the tissue (skin and fat) thickness varied at P_elect = 120 W. The peak pressures exhibit a nearly constant trend when the thickness of skin or fat tissue is altered, except for the case of skin tissue at f > 1.1 MHz. These results suggest that the thickness of skin and fat have little influence on the pressure in the blood region.

Figure 6.

Peak pressure changes with the tissue thickness. (a) and (b) correspond to change of skin and fat layers, respectively.