Prediction of HIFU Propagation in a Dispersive Medium via Khokhlov–Zabolotskaya–Kuznetsov Model Combined with a Fractional Order Derivative

: High intensity focused ultrasound (HIFU) has been proven to be promising in non-invasive therapies, in which precise prediction of the focused ultrasound ﬁeld is crucial for its accurate and safe application. Although the (KZK) equation has been widely used in the calculation of the nonlinear acoustic ﬁeld of HIFU, some deviations still exist when it comes to dispersive medium. This problem also exists as an obstacle to the Westervelt model and the Spherical Beam Equation. Considering that the KZK equation is the most prevalent model in HIFU applications due to its accurate and simple simulation algorithms, there is an urgent need to improve its performance in dispersive medium. In this work, a modiﬁed KZK (mKZK) equation derived from a fractional order derivative is proposed to calculate the nonlinear acoustic ﬁeld in a dispersive medium. By correcting the power index in the attenuation term, this model is capable of providing improved prediction accuracy, especially in the axial position of the focal area. Simulation results using the obtained model were further compared with the experimental results from a gel phantom. Good agreements were found, indicating the applicability of the proposed model. The ﬁndings of this work will be helpful in making more accurate treatment plans for HIFU therapies, as well as facilitating the application of ultrasound in acoustic hyperthermia therapy.


Introduction
Although pioneering clinical studies of focused ultrasound were carried out as early as the 1940s [1,2], it did not attract intensive research interest until the end of the 20th century and the beginning of 21st century, during which several theoretical models were developed, improved and then broadly accepted [3][4][5][6][7][8][9]. In the past decades, high intensity focused ultrasound (HIFU) has played an increasingly significant role in the study of non-invasive therapies by demonstrating unique advantages in safety, effectiveness and high efficiency [10][11][12][13][14]. However, the applications of HIFU are still limited, and clinical treatments are only available for limited sites [11,[14][15][16].
One of the challenges confronting HIFU treatments is the spatial precision of tissue ablation. Several techniques such as ultrasound B-Scan and Magnetic Resonance Imaging (MRI) have been combined with HIFU to achieve real-time monitoring of focal areas [17][18][19][20]. With these methods, the actual focal profiles were usually found to deviate from those predicted through theoretical models [21,22]. As was indicated by Petrusca et al., the shift of the focal point away from the prescribed α ∝ ω y , y = 1 − 2. For most tissues, the attenuation factor y sits in the range from 1 to 1.7 [32,33]. However, in the KZK model, the propagation loss caused by viscosity and thermal conduction is considered to be proportional to the square of the working frequency (y = 2), which is actually only valid for fresh water [32,33]. Furthermore, y = 2 for fresh water causes a third derivative term to appear in the KZK equation. The third derivative as well as the derivative of acceleration have not yet been well clarified in the physical overview. In the context of Newton's law of motion, acceleration is directly affected by the force, and the derivative of acceleration has no obvious physical meaning. Secondly, sound velocity, which is treated as a constant in the KZK equation, usually changes with frequency in biological tissues. In HIFU fields, propagation nonlinearity could be non-negligible due to the high acoustic pressure, and the orders of harmonic waves could be rather high in many cases [5][6][7]9,29]. Therefore, the descriptions of both attenuation and sound velocity should be modified when biological tissues are present on the wave path.
To overcome these shortcomings in the existing models, efforts need to be made to account for the case of a non-integer power index in the Kramers-Kronig relationship. In the frequency domain, with the help of the classic Laplace Transform, the acoustic field could be easily accessible, but only valid for linear models. The Rayleigh model proposed by Wojcik et al. in 1995 was also not applicable in cases using wide-bandwidth acoustic pulses [34]. The fractional order derivative method described by Makris and Constantinou [35] was believed to be an effective approach to solve this problem. However, the complexity in mathematics has made numerical analysis hard to achieve, hence restrained its application in developing HIFU theories. After that, Szabo et al. utilized a convolution integral and further developed the theory of the fractional order derivative, in which Fourier Transform was adopted to define a derivative of non-integer order, while the basic principles of the Kramers-Kronig relationship was still preserved [36][37][38]. Some other calculation algorithms have also been proposed. Treeby et al. developed a popular, open source nonlinear simulation tool named the k-wave toolbox to simulate nonlinear wave propagation [39]. Prieur et al. proposed time-fractional acoustic wave equations [40]. Inspired by these algorithms, simulation of the KZK equation has made definite progress in the past few years.
In this paper, a fractional order derivative was introduced to modify the KZK model to address the power-law relationship (non-integer power index), and a frequency-dependent sound velocity was employed to account for the dispersive behaviors of media. Numerical and experimental verifications were carried out to demonstrate the different behaviors between the modified model and the original KZK model. Results showed that the obtained modified KZK (mKZK) equation is in better agreement with the experimental results. The findings in this paper will promote not only the prediction and design of focused sound fields and acoustic transducers, but also the therapeutic applications of ultrasonic treatments.

The KZK Equation
The KZK equation is an extended form of Burgers model, and an approximation of the Westervelt equation, written as [41] ∂ 2 p ∂z∂τ in which p is the acoustic pressure, c 0 is the sound velocity, δ is the sound diffusivity, β is the nonlinearity coefficient, ρ 0 is the ambient density of the medium, ∆ ⊥ is the transverse Laplace operator (defined as ∂ 2 ∂θ 2 in cylindrical coordinates), and τ = t − z/c 0 is the time delay at the axial distance of z with t being the time. The three terms on the right side of Equation (1) represent the diffraction, attenuation and nonlinearity, respectively. The attenuation caused by viscosity and thermal conduction is hence proportional to the square of the angular frequency ω, i.e., α(ω) = ω 2 δ/(2c 3 0 ). However, in biological tissues, the attenuation factor is usually a non-integer less than 2, the attenuation term is hence only valid for describing wave propagation in fresh water (attenuation factor y = 2). If the attenuation parameter is directly set for fresh water rather than considering the actual properties of the media, the accuracy of calculation would certainly be undermined. Furthermore, the sound velocity c 0 here is regarded as a constant for all harmonic components, which is in conflict with the inherent Kramers-Kronig dispersion relations where the phase velocity varies with increasing frequency.

The Modified KZK Model
By introducing Fourier Transform, the nth time derivative of time-dependent acoustic pressure p(t) could be written as Here F + and F − represents operators of the Fourier transform and its inversion, respectively. As the order number y is a non-integer, the fractional order derivative is then defined as a convolution [38], where Γ(·) is the Gamma Function, i.e., Γ(x) = ∞ 0 ξ x−1 e −ξ dξ. The definition in Equation (3) is hence the Riemann-Liouville fractional derivative [42]. Therefore, the fractional order derivative is not only determined by the pressure value at the time point t, but is also related to its history range from −∞ to t.
Since the attenuation factor y is a non-integer in lossy media, the wavenumber k could then be considered as a complex, while its squared form in the low-frequency approximation is in which i is the imaginary unit, and (−i) y can be expressed with trigonometric functions as (−i) y =cos(yπ/2) − isin(yπ/2). With the wave number k = ω/c 0 + iα 0 (−iω) y /[cos(yπ/2)], the phase velocity can be calculated as where Re(·) represents the real part of the complex value.
Considering the non-integer attenuation factor y and the dispersion of phase velocity, the classical KZK equation could be modified as where the attenuation term and the phase velocity c(ω) could be calculated according to Equations (4) and (5), respectively. Consistent conclusions can be found between Equation (6) and the work of Zhao et al. [43], which is an extension of an earlier model for ultrasound propagation in power-law media proposed by Kelly et al. [44]. In addition, Equation (6) will not hold if y = 1. However, in the framework of discussing HIFU propagation problems, the y = 1 case is of no importance [37,38] and is not of interest here.

The Numerical Algorithm
In the numerical analysis, both the KZK model and its modified form were solved with the finite difference time domain (FDTD) method for the acoustic field emitted from a single-element self-focusing transducer. For the traditional model, coordinate transformation Z = z/F, R = r/a, T = ωt, P = p/P 0 were introduced, with F, r, a, P 0 being the geometrical focal length, the radial coordinate, the aperture radius of the transducer, and the surficial acoustic pressure, respectively. The following assumptions were then made, in order that the KZK equation could be normalized to the following form, The values of the physical constants used for acoustic modeling were ρ 0 = 1000 kg/m 3 , c 0 = 1486 m/s, β = 3.5, α = 0.025 Np/m at 1 MHz, and µ = 2 for fresh water [45]. The FDTD algorithm adopted here was generally the same as that used in a previous study [45,46], in which Equation (10) was decomposed into three independent equations, accounting for the diffraction, attenuation and nonlinearity, respectively. Based on an orthogonal spatial grid, the discretized forms of these equations were expressed as, And Here i and j were the spatial coordinate indexes in the radial and axial directions, respectively, and n was the time step, i.e., P n i,j = P(i · dr, j · dz; n · dt). In the simulations for the mKZK equation, the major difference was the differential form of the fractional derivative in Equation (6), where The FDTD simulations were then accomplished through a self-developed FORTRAN code package running on a ×64 PC platform. In the calculation, the boundary condition was symmetrical and had equal amplitude acoustic pressure driving conditions, which was also the boundary condition generally used to calculate HIFU [29,41].

Experimental Setup
The thickness of each phantom sample was L s = 42.3 mm, with a density of ρ 1 = 1000 kg/m 3 so that it could stably suspend in the water. The distance between transducer and the phantom was 37.7 mm. Following the same protocol used in our previous work [45], the nonlinearity parameter was measured as β = 4.2. Since the only difference in gel recipe between the two works is the introduction of amino polystyrene microspheres in this paper, the same β value indicates that microspheres did not influence the nonlinear propagation of waves. To confirm this, we measured the ratio between the second harmonic and fundamental components, and found it was identical to that in [45] under the same sonication conditions. The attenuation coefficient and sound velocity of phantom, as functions of frequency, were measured through a broadband spectrum method [48,49]. In the measurement, two planar piston transducers (Immersion, Unfocused, Panametrics, Waltham, MA, USA) calibrated with a needle hydrophone (HNC-1000, ONDA Corp., Sunnyvale, CA, USA), were placed on the opposite sides of the phantom, with one of them driven by a broadband pulse generator (5900PR, Panametrics). The reflected and transmitted acoustic signals were then acquired by the transducers and digitalized with a digital oscilloscope (54830B, Agilent, Santa Clara, CA, USA). (Device connection was similar to that in [49]). In the measurements, 8 continuous pulse sequences were acquired and averaged to reduce the signal-to-noise ratio (SNR).
Prior to the HIFU experiments, a low-level driving voltage was used to drive a customized HIFU transducer (Chongqing Haifu Med. Tech. Co., Ltd., Chongqing, China) to emit a linear sound field. Simultaneously, effective parameters of the transducer, such as effective radius, radiation profile, and angle of divergence were obtained by adjusting the transducer parameters in the KZK calculations, such that the linear field predicted via KZK was consistent with that measured [45]. The effective parameters were then used in the simulations of both the KZK and mKZK models. As a result, the HIFU transducer (working frequency 1.12 MHz) had an effective aperture radius of 48.6 mm and a geometrical focal length of 101.5 mm. As illustrated in Figure 1, the transducer was immersed in water and driven with signals from a signal generator (33250A, Agilent, Santa Clara, CA, USA) amplified by a broadband power amplifier (2200L, E&I, Rochester, NY, USA). The input voltage was set as 465 mV (20 cycles; burst period, 10 ms; duty cycle, 0.18%). Another needle hydrophone (HNA-0400, ONDA Corp., Sunnyvale, CA, USA) was mounted on a customized three-dimensional (3D) scanning system (Controller Model: XPS-C8, Newport, CA, USA) to scan the HIFU field. To suppress possible acoustic cavitation in surrounding liquid, the water was processed with a self-developed water degassing and deionizing system. The temporal and spatial scanning procedure was controlled via the GPIB interface (National Instruments, Austin, TX, USA).

Non-Dispersive Water
In order to verify the validity of the modified model, the mKZK equation was used to predict the sound field distributions generated from the transducer. In this case, the media in the direction of propagation was degassed water, and the surface pressure of the transducer was set to be 0.4 MPa. Then, the results were compared with those obtained from experimental measurements as well as from the original KZK model. The results are presented in Figure 2 for comparison. Due to the axial-symmetry of the sound field, only the sound field distribution in the axial direction was studied. Also, since the major concern of this paper is to investigate how the focal-shift could be accurately predicted, the pressure profiles are all presented in a normalized way, so that the focal-shift effect is more intuitive and easy to observe. For the total pressure distributions presented in Figure 2a, the axial distribution of acoustic pressure seems identical for the modified and original KZK models, which provides reliable proof that the current theoretical modification did not compromise the accuracy of the sound field prediction in non-dispersive media. With the help of fast Fourier transform (FFT) algorithm, further analysis was then carried out by decomposing the total sound pressure into the superposition of linear and nonlinear components. In Figure 2b-d, all the components exhibited good agreement between the results calculated from the modified and the original KZK models, showing that both models are applicable for sound field prediction under the experimental conditions mentioned above. It should also be noted that, in Figure 2 the locations of the pressure peaks from the two models were exactly the same for all components, although the measured axial beam-widths seem a bit narrower than both theoretical predictions, especially for harmonic components. Since the linear fields were found to be almost identical for the measured and predicted results, we speculate that some far-field attenuation factors such as dissolved oxygen in water, or bubbles might exist. However, this does not affect the conclusion on the location of the maximum pressure.

Dispersive Phantom
To incorporate the phantom model into the study, the acoustic parameters of the phantom was firstly characterized according to procedures described in earlier studies [48,49]. In brief, the acoustic phase velocity c(f ) inside the phantom material was determined by [48] where the sound velocity in water c w was considered as 1500 m/s. When the sound velocity inside the phantom was measured, the acoustic signals p s and p f were the acoustic pressure acquired by a transducer before and after the phantom was inserted into the acoustic path. While p 1 and p 2 were the pressure of the reflected signals from the first and second water/phantom interfaces, respectively, θ w (f ), θ s (f ), θ 1 (f ) and θ 2 (f ) were the corresponding phase spectra. The frequency dependence of the attenuation coefficient α(f ) was calculated according to [48,49] where A w , A s , A 1 and A 2 were the amplitude spectra corresponding to the above-mentioned phase responses. Figure 3 plots the measured attenuation coefficient and acoustic velocity as a function of frequency. In the frequency range 1.5-3.1 MHz, the sound velocity increased by about 26 m/s. The acoustic attenuation coefficient increased from 0.59 dB/cm at 1.5 MHz to 1.79 dB/cm at 3.1 MHz, indicating the attenuation factor being y = 1.83 for the phantom. The attenuation factor was obtained from the exponential fitting using a curve fitting toolbox in Matlab. It should be noted that sound velocity could fluctuate due to the thermal effect of focused ultrasound. However, in this work the duty cycle was as low as 0.18%, and a medium-level surface pressure of up to 0.4 MPa was chosen for the transducer. Thus, no significant temperature elevation was observed during the experiments, and the thermal-induced change in sound velocity could be neglected [50]. It should also be mentioned, that to account for the thermal-effect, an appropriate bio-heat transferring equation should be incorporated with the current model. The main concern then lies in the dispersion due to microsphere scattering. In clinical applications, the major dispersion originates from tissues like fat, whose attenuation coefficient is generally larger than that of body tissue. Therefore, in comparison with other research in which the parameters of the phantom were nearly the same as body tissue (e.g., liver and spleen), the choice of phantoms with larger attenuation coefficients in the present experiments might provide results in better agreement with the actual situation. Meanwhile, since the measured physical parameters are close to those of dense fat, this setup could be regarded as a simplified mimic of the abdomen in HIFU therapies. As is illustrated in Figure 4, the acoustic field distribution along the axial direction is examined by sitting the phantom on the acoustic path of the transducer. The results show the comparison between measured data and simulated results obtained from both the modified and traditional KZK models. Figure 4a describes the total pressure and Figure 4b-d shows the fundamental, second harmonic and third harmonic components, respectively. For the total pressure distribution, it is clearly observed that the peak-pressure location predicted through the mKZK model agrees well with that acquired in experiments, while the data calculated from the traditional KZK model show a deviation from the previous two groups. Note that the axial peak-pressure deviation is quite small for the fundamental components, but it gradually becomes significant when more harmonic components appear. Thus, the overall peak-pressure deviation observed in Figure 4a is mainly caused by higher-order harmonic components, indicating that the significance of the theoretical modification relies on high nonlinearity, such as in HIFU. This phenomenon could be addressed with the results shown in Figure 3, where higher-order harmonic waves that occupy higher frequency bands exhibit more notable sound velocity/attenuation dispersions, thus play a more dominant role in the modification of the KZK model. Meanwhile, it can be seen from Figure 4 that the modified equation has larger variation than the KZK equation. The effect of the mKZK equation on the simulation results can be thus clarified as providing more precise prediction for experimental results. The results in Figure 4 give persuasive proof that theoretical modifications made previously are necessary and valid. The beam narrowing effect caused by data normalization still exists. However, in actual treatment more attention is paid to the location of the focus point, because the location determines the heat distribution area, which significantly alters the biological properties of the treatment area. Among existing studies, most researchers have focused on the thermos-lensing effect [23] and bubble formation induced focal region distortion [25,26], while some also mentioned the acoustic radiation force induced tissue displacement [27]. With the results presenting the dispersion of sound velocity and attenuation in tissues, the deviation of focal spots can be well explained in combination with the strong nonlinearity nature of HIFU.

Dispersion-Induced Focus Shift
It is of great importance to evaluate how the mKZK model demonstrates its significance in HIFU applications. In Figure 5, comparisons are carried out to display how the wave distribution behaves differently before and after inserting the phantom sample into the wave propagation path. It can be observed that, for either data from experiments or from the mKZK model, although only a slight difference is seen in the axial wave profile when examining the fundamental components, an axial focus shift is more evident for the overall acoustic pressure since it is highly affected by the harmonic components. It is explained here that, due to the dispersive nature of the phantom, higher-order harmonic components require larger values of both sound velocity and attenuation coefficient in the KZK model, urging the overall wave profile to move forward to the transducer. In the present study, the focus shift distance was quantified for the overall axial acoustic pressure distribution and its decomposed components, and this is listed in Table 1. The focus shift, which is found to be of millimeter magnitude for the studied condition, cannot be ignored, especially for clinical HIFU studies. On the one hand, in typical HIFU applications, the much higher surface pressure would induce even stronger acoustic nonlinearity in the focus area, and up to tens of orders of harmonic components might contribute to the overall acoustic responses. In that case, dispersion in acoustic velocity as well as the attenuation coefficient could cause larger focus shifts. On the other hand, even for the millimeter-level focus shift exhibited with the current setup, an impressive amount of acoustic energy would be deposited outside the designated focus area. As demonstrated in Figure 6, about 25% of the −3 dB focal region would fall outside of the one predicted by the traditional KZK equation, when the phantom sample is introduced into the transducer axis. In addition, the absence of shockwaves in the experiment should be noted, and this indicates that the influence of shockwaves could be ignored in the theoretical model under such experimental parameters.

Discussion
In previous studies, shifts between the predicted and observed focal regions have been frequently reported in HIFU-related studies. Usually, the shift was measured to be around 1-3 mm for 1-MHz HIFU transducers [24], 4-5 mm for~2.2-MHz excitation [51], and an empirical formula was also used to calculate the focus shift [22]. Although different possible mechanisms have been proposed and a series of correlation studies have been carried out [21][22][23][24][25][26][27], detailed theoretical proof that could quantitatively explain the inherent mechanisms of the observed focal shifts is still lacking. The proposed model modification here clarifies how the acoustic dispersion played a role in the complicated physics of this problem.
The experiments carried out here have demonstrated that the observed focal shift should come from the dispersive behavior of sound velocity and attenuation in media. For different harmonic components in the HIFU beam, their sound velocities could be different at the gel/water interface. Speculating from Snell's law, these different components could actually propagate along slightly different paths in the phantom, causing the focal point to shift its location. That is also why the main difference between the two models is observed for the harmonics in Figures 4 and 5-because the sound velocity of the harmonics deviated further away as their frequencies were higher. During this process, the attenuation should also have contributed in a dispersive way. However, the influence of dispersive attenuation could not be separated from that of sound velocity.
Despite the general recognition that the Westervelt equation and SBE provide higher precision in non-dispersive sound field prediction than KZK, difficulties in predicting accurate sound field distribution in strong dispersion media still remains a great challenge [22]. Moreover, due to the existing complexity of the Westervelt equation and SBE, adding more modifications to these two equations could be over-whelming and/or time-consuming for clinical applications that require real-time monitoring. Considering that KZK shows a balance between accuracy and computation burden, this paper chose to further modify the KZK equation as a straightforward way to achieve improved simulation of HIFU propagation, especially to work out how the focal-shift could be predicted accurately. In this work, for consistency of the experiment, the mKZK and KZK in ordinary media was first confirmed. Then, when measuring the parameters in the strong dispersion medium, mKZK gave a more accurate prediction of the focal shift. Thus, by using the mKZK equation we can quickly and accurately predict the HIFU field in complex media.
However, further theoretical studies are still needed since the current model is unable to eradicate all the possible unfavorable factors in predicting the characteristics of HIFU. For example, defocusing effects or shifts in focal position might be caused by the layered tissue effect, where sound speed/attenuation may vary in different tissues layers. The thermal lesion effect, in which the tissue properties change due to the heating of HIFU, could be more difficult to include in the modeling. A full solution could be even more challenging if other possible mechanisms, including acoustic radiation force, acoustic cavitation and inconsistent thermal deposition [27] are also considered.
The focal shifts observed in phased-array-based HIFU devices are also notable. For instance, a focal shift of about 2-mm was observed along the transducer axis in multiple-layered soft tissues sonicated with a 65-element phased array transducer [22]. To overcome this problem in phased-array HIFU, possible solutions could be obtained by drawing lessons from the underlying mechanisms discussed above.

Conclusions
Although HIFU technology based on the KZK equation calculation has been widely accepted and used in the clinical setting and transducer designs, the absence of an accurate theory to predict the sound field inevitably limits the application of the ultrasound focusing. In this work, a mKZK equation is proposed to predict the HIFU field established with a spherical focusing transducer. Meanwhile, the accuracy of the methods is verified through experiment. This method could improve the computational accuracy of the KZK equation in dispersive media, which is similar to human tissue. Simulation and experimental results show that the focus area will shift towards the transducer and the offset increases as the nonlinearity becomes higher. Therefore, in the process of HIFU transducer design, the impact of dispersion on the results need to be taken into account, in order that accurate sonication can be achieved. By modifying the KZK equation, the findings will help with transducer design and the application of the HIFU. This will also help to ensure the stability and safety of HIFU and further accelerate its clinical applications.