Time-domain simulation of ultrasound propagation with fractional Laplacians for lossy-medium biological tissues with complicated geometries

Simulations of ultrasound wave propagation inside biological tissues have a wide range of practical applications. In previous studies, wave propagation equations in lossy biological media are solved either with convolutions, which consume a large amount of memory, or with pseudo-spectral methods, which cannot handle complicated geometries effectively. The approach described in the paper employed a fractional central difference method (FCD), combined with the immersed boundary (IB) method for the finite-difference, time-domain simulation. The FCD method can solve the fractional Laplace terms in Chen and Holm's lossy-medium equations directly in the physical domain without integral transforms. It also works naturally with the IB method, which enables a simple Cartesian-type grid mesh to be used to solve problems with complicated geometries. The numerical results agree very well with the analytical solutions for frequency power-law attenuation lossy media.


I. INTRODUCTION
Ultrasound has been widely used as a diagnostic tool in medical imaging for over 50 years.Its applications have been extended to include destruction of kidney stones, killing malignant tissues, and cosmetic surgery (Duryea et al., 2014;Maxwell et al., 2015).Focused ultrasound beams can also be used to remove the brain tumors (Hill and ter Haar, 1995).A recent study shows that ultrasound techniques in conjunction with microbubbles can be used to safely open the blood-brain barrier (BBB) for brain drug delivery, which is a new treatment in instance of stroke or Alzheimer's disease (Konofagou, 2012;Wu et al., 2015).Therefore, using ultrasound safely requires accurate planning, which motivates accurate ultrasound simulation techniques.
Due to the effects caused by heterogeneous tissue media, thermal conduction, viscous dissipation, and chemical relaxation processes, ultrasound propagation processes inside human tissues or lossy media are much more complicated than those in simple media such as air.A lossy medium is a medium in which a significant amount of acoustic energy is absorbed per unit distance traveled by a sound wave.Most biological tissues can be considered lossy media.Sound attenuation is usually used to quantify an energy loss in lossy media, which has also been found following the power law in the frequency domain as (Blackstock, 1985;Szabo, 1994) a ¼ a 0 jxj y ; (1) where a is sound attenuation with a unit of Np/m, a 0 is absorption coefficient, x is angular frequency, and y is power law exponent of the specific material with a value between 0 to 2.
To simulate sound propagation inside lossy media, different numerical approaches have been proposed over the past few decades.Classical thermo-viscous theory predicted that the acoustic sound would attenuate inside lossy media due to thermal conduction effects and viscous dissipation (Morse and Ingard, 1968).But the theory can only predict square-law attenuation (y ¼ 2).Sparrow and Raspet implemented finite difference time domain method (FDTD) in 2D axisymmetric domain with nonlinear effects where the power attenuation law was not considered (Sparrow and Raspet, 1991).Szabo (1994) proposed a time domain causal convolution operator that account for both power law absorption and dispersion ð1 < y < 2Þ.However, the implementation of convolution integration consumes extremely large memory (Norton and Novarini, 2003), especially when the computational domain is large.To improve the efficiency of the calculation, a fractional Laplacian model replaced the convolution model (Chen and Holm, 2003).Treeby and Cox implemented this model with the pseudo-spectral method, but theirs "k-space" method still needs to transform spatial derivatives back to the frequency domain, which makes it difficult to deal with complicated geometries (Treeby and Cox, 2010).
This paper presents an approach in which a simple structured Cartesian grid mesh can be used to solve ultrasound propagation problem with any irregular geometry of lossy media describable by the frequency power attenuation law, Eq. ( 1).The finite-difference time-domain (FDTD) method is coupled with the IB method to accommodate complicated geometries (Xu et al., 2011;Ke and Zheng, 2015;Ke and Zheng, 2016).The lossy medium is modeled with the Chen a) Electronic mail: jjzhang@ku.edu and Holm's equation (Chen and Holm, 2003).In order to calculate the fractional Laplacian terms in the model, the fractional central difference method (FCD) (Liu et al., 2015) is used.The perfectly matched layer (PML) boundary is used to mimic a free space condition (Hu, 1996;2005).This new approach, different from the "k-space" method, does not need additional correction factors or integral transforms and can accommodate complicated geometries with a simple structured mesh.
The governing equations will be presented in Sec.II, along with the detailed numerical method.Simulation examples of sound propagation inside a lossy medium are given in Sec.III.After comparing the numerical simulation results with the analytical power attenuation law, we will discuss the results and offer the conclusion in Sec.IV.The order of accuracy of the scheme is also evaluated and presented in the Appendix.

II. NUMERICAL METHOD A. Governing equations in lossy media
Two media are considered in this study: water and the lossy-medium biological tissue.The linearized Euler equations for wave propagation in water are where u 0 ; p 0 ; and q 0 are time averaged velocity, pressure, and density of water.u, p, and q are their corresponding acoustic fluctuations, c 0 is the speed of sound in water, and f p is the fictitious term for implementing the immersedboundary (IB) method, which will be explained later.In water, if there is no background flow and the background pressure is a constant value, Eqs.
(2) and (3) become In lossy media, conservation of mass is expressed as where q 1 is the density of the lossy medium.Conservation of momentum is still in the form of Eq. ( 4).The relation between acoustic pressure and acoustic density in the lossy medium is derived by Treeby and Cox (2010) based on Chen and Holm's equation as q; (7) where c 1 is the speed of sound in the lossy medium, and s and g are proportionality coefficients.Substituting Eq. ( 6) into Eq.( 7) gives the acoustic pressure propagation equation in the lossy medium as The proportionality coefficients, s and g, follow the relations where the coefficients are chosen as y ¼ 1:9, a 0 ¼ 2:9858 Â10 À10 dB Hz Ày m À1 , c 1 ¼ 2000 m/s, q 1 ¼ 1500 kg/m 3 to simulate a relatively hard tissue.The corresponding proportionality coefficients based on Eqs. ( 9) and ( 10) are

Computation of fractional derivatives
Consider the Riesz fractional derivatives equal to the fractional Laplacian operator (Yang et al., 2010) given below, The fractional central difference (FCD) method (Liu et al., 2015;Qin et al., 2018) can be used to numerically approximate the Riesz fractional derivatives in Eq. ( 11).For a continuous function f ðx; zÞ, each term in Eq. ( 11) can be expressed as where the weight function x q is The fractional derivative terms in Eq. ( 8) are approximated with the summation of the weighted pressure or velocity gradients in space using Eq. ( 12).Comparing this method with the convolution method used in Norton and Novarini (2003), the memory requirement is reduced significantly.
For an irregularly shaped lossy medium, the number of computational grid cells are different in the x-and z-directions.Therefore, it is necessary to identify the range of the lossy medium in each direction.A horizontal line crosses the boundaries in Fig. 1 can have two intersections at the left and right boundaries of the geometry.Those intersection points are not necessarily located on the boundaries because of the Cartesian grid mesh used in the simulation.Therefore, they are approximated with the closest grids, eðjÞ and wðjÞ, near the boundaries.Similar approximation is needed along a vertical line as shown in Fig. 1.When the IB method is implemented, those boundary grids, eðjÞ; wðjÞ; nðiÞ; and sðiÞ, are flagged and stored at the beginning of the computation, as demonstrated in Fig. 1.This makes the implementation of FCD method very simple and efficient.

C. Immersed-boundary method
The f p in Eq. ( 3) is the fictitious term in the immersedboundary method, which is used to represent the material change between water and the lossy medium.The existence of the boundaries of the lossy medium is represented by switching the fictitious force in the equation.The computation can be performed as if there are no boundaries.In Eq.
(3), the fictitious force term for the immersed-boundary method implementation is expressed as outside the lossy meidum; ðr Á uÞ; inside the lossy medium: The acoustic pressure outside the lossy medium is modeled with original pressure equation, Eq. ( 5), as the linearized Euler equation with f p ¼ 0. The convection terms are neglected assuming the background flow speed is low.The pressure inside the lossy medium is modeled with Treeby and Cox's equation as ðr Á uÞ.This is constructed so that when substituting f p back to Eq. ( 3), the lossy medium model equation, Eq. ( 8), is resumed.Therefore, by switching the fictitious force term, the same computational solver can be used simultaneously for both materials.Moreover, with this method, a Cartesian grid mesh, regardless of the complicated object boundaries in the simulation domain, can be used to solve for the acoustic field.
To absorb numerical reflections by computational domain boundaries, the PML boundary (Hu, 1996;2005) methods are used on the outside boundaries of the domain.The PML boundary condition has been successfully implemented and verified in our previous work (Zheng and Li, 2008;Ke and Zheng, 2015;Ke and Zheng, 2016) for the FDTD simulation.

III. NUMERICAL SIMULATION AND RESULTS DISCUSSION
A. Ultrasound propagation in the lossy medium and water Our first simulation is for acoustic propagation in the lossy medium only.A simulation for ultrasound propagation in water is also conducted as a reference.The size of the computational domain is 0.06 m Â 0.06 m.A uniform Cartesian-type mesh with the grid size of Dx ¼ Dz ¼ 2:5 Â10 À5 m ensures at least 20 grid points are used per wavelength for ultrasound frequencies up to 1.5 MHz.The speeds of sound of water and the lossy medium are set to 1500 and 2000 m/s, respectively.The density of water and lossy medium are set to 1000 and 1500 kg/m 3 .The Courant-Friedrichs-Lewy (CFL) number is chosen to be 0.3 to satisfy a stable computation.It should be noted that if we reorganize Eq. ( 8) by moving the convection term to the left hand side, the left hand side is still in the form of Euler equation.The two fractional Laplacian terms on the right hand side of equation are calculated as the summation of a series, which is not part of the finite difference scheme.Therefore, the original stability criteria of the scheme in Zheng and Li (2008)  right, top, and bottom boundaries of the computational domain, respectively.Four receivers are placed along the centerline of z ¼ 0:03 m, at x ¼ 0:005, 0.015, 0.025, and 0.035 m to record the acoustic pressure histories in these locations.It should also be noted that we used this 1D problem to perform 2D simulation, for the purpose to verify the implementation of PMLs in the z-direction for preparing for the following IB method study.It is evident that the implementation of PMLs was successful.More details on the PML implementation can be found in (Zhang et al., 2016).
Figure 2 shows acoustic pressure contours at simulation time t ¼ 22:5 ls.The plane wave in lossy medium travels further than that in water at this moment.However, the pressure magnitude in water is much higher at this moment, which can also be observed from receiver pressure histories in Fig. 3.The phase shift in Fig. 3 between solid and dashed lines can be explained by the speed of sound difference in the two media.Pressure histories of the four receiver locations in water show almost the same magnitude, while the pressure magnitude is gradually decreasing along the wave propagation direction in the lossy medium.The phenomenon agrees with the definition of lossy medium that the acoustic energy is absorbed with the distance travelled in the medium.

B. Comparison with the analytical solution
To quantitatively verify the simulation results, a comparison with the frequency-domain power attenuation law is conducted.The analytical solution of attenuation in the lossy medium is described in Eq. (1).To obtain the attenuation from the time-domain numerical calculation, the pressure histories recorded in Sec.III A are used.After applying the Fourier transform, the power spectrum density (PSD) can be used to calculate the sound attenuation (SA) in the frequency domain as By averaging the sound attenuation over the wave propagation distance d, the normalized sound attenuation a from the numerical simulation is calculated as It should be noted that the unit of a in Eq. ( 1) is Np/m.A unit conversion from Np to dB is necessary to compare with the SA obtained from Eq. ( 16).The comparison of sound attenuation values between the numerical solutions and the analytical solutions for different power law exponents, y, are plotted in Fig. 4. All simulation results show near-perfect agreement with the analytical solutions.

C. Simulation of ultrasound propagation in a complicated geometry
The successful verification of the FCD in the lossy medium simulation enables further implementation for simulating wave propagation with more complicated geometries.A piece of ring-shaped lossy medium, which is intended to represent a bone-type biological material, is placed in the middle of the computational domain.The inner diameter of the ring is 0.05 m and the outer diameter is 0.015 m.The ultrasound source is placed on the left boundary from z ¼ 0.025 to 0.035 m to represent an ultrasound transducer.The total simulation time is 38 ls, which allows the incident ultrasound waves to fully pass through the ring-shaped lossy medium.
The detailed simulation setup is shown in Fig. 5.The areas in blue, grey, and white indicate water, lossy medium, and PML, respectively.The simulation domain is surrounded by PMLs with a thickness of 0.01 m to eliminate reflection waves from the computational boundaries.The dashed line on the left is to indicate that the area on the left side of the line will be switched into a PML once the wave has fully passed through.The medium interfaces, represented by the two concentric circles in the domain, are not aligned with the Cartesian grid mesh.Therefore, the immersed-boundary method is used to accommodate the cut-through Cartesian grid mesh between the interfaces of water and the lossy medium.The closest node points are used as the boundary points in the case of cut-through.
Figure 6 presents acoustic pressure contours recorded at four different time moments.Reflected waves can always be found at the interface, which are formed by the material difference of the two media and the geometrical shapes of the interfaces.For example, when waves propagate from the lossy medium to water in Fig. 6(b) and from water to lossy medium in Fig. 6(c), the two plots clearly capture the backward wave reflections.When the waves leave the lossy medium in Fig. 6(d), they are much weaker than the incident waves in Fig. 6(a), which shows the combined effects of wave reflection and dissipation due to the lossy medium.
To illustrate the capability of the simulation, a multiplesource simulation is conducted.The same type of plane waves are specified on top, bottom, left, and right boundaries.Pressure contours at different simulation moments can be found in Fig. 7.The wave patterns are more complicated due to interactions of multiple incoming waves with the ring-shaped lossy medium.Comparing Fig. 6(c) with Fig. 7(d), which are at the same moment, the effect of multiple incident sources can be observed clearly.

IV. CONCLUSION
When sinusoidal waves propagate into a lossy medium, they experience dissipation caused by viscosity, heat conductivity, and relaxation process.Therefore, incident acoustic waves are attenuated by those processes, which results in a power-law attenuation in the frequency domain.This phenomenon can be modeled and numerically simulated with the FDTD method.When solving for fractional Laplacian derivatives in the lossy medium acoustic propagation equations, the FCD method is employed (Liu et al., 2015;Qin et al., 2018).The procedure, which is based on a relatively simple algorithm does not require large memory involve integral transform.The implementation of the FCD combining with the IB method maintains the second-order spatial accuracy.
The numerical simulation results show a gradually reducing pressure level along the wave propagation direction in the lossy medium and agree with the analytical solution of a power attenuation law.The implementation of the immersedboundary method enables ultrasound propagation around irregular geometries to be simulated with a simple structured Cartesian mesh.Moreover, this method can be used to solve multi-media, multi-source problems.Since the mesh is already designed to resolve very high frequency signals (over 1 MHz) for ultrasound propagation problems, the density of the mesh usually enables accurate interpolation near the boundary of the medium geometry unless there are extremely small local features.For very small local geometrical features, we can increase the local mesh resolution to increase the accuracy of the immersed-boundary method.This capability for handling complex geometries with multiple material interfaces overcomes the limitation of the previous work in the area of lossy medium simulation.Furthermore, with the same approach, the method can be readily extended to simulate 3D problems.

ACKNOWLEDGMENTS
This research is partly supported by the U.S. Army under a cooperative Agreement No. W911NF-14-2-0077.

APPENDIX: NUMERICAL SCHEME ORDER OF ACCURACY VALIDATION
While a second-order accurate scheme is used in the simulation for both time and space (Zheng and Li, 2008), the approximation of boundary treatments used in the FCD and IB methods can possibly introduce additional errors in spatial discretization.Therefore, it is necessary to re-evaluate the order of spatial accuracy after implementation of the IB and FCD.
Three cases are tested here, with the water-only simulation as a baseline.The lossy-medium simulation is conducted to investigate the accuracy of the FCD method.Then, the simulation with the ring-shaped lossy medium is used to investigate the accuracy of the implementation of the IB method combining with the FCD method.Four levels of grid sizes are used in each simulation case, namely, coarse, medium, fine, and finest.The grid size h of each level is 5 Â 10 À5 , 2:5 Â 10 À5 , 1:25 Â 10 À5 ; and 0:625 Â 10 À5 m, respectively.A very small time-step, Dt ¼ 1:25 Â 10 À9 s, is used to ensure the stability of the simulation.The simulation results at the moment of t ¼ 7:5 ls, which is when the wave front has fully passed the center of the domain, are used.All the nodal point pressure values along the centerline z ¼ 0:015 m are used to calculate the L 2 -norm of the acoustic pressure.The L 2 -norm error is calculated as Since the exact solution is not available, the finest grid solution at h ¼ 0:625 Â 10 À5 m is used as the reference.Table I gives the L 2 -norm errors of the different grid resolutions in the three cases.The observed order is computed using P ¼ logðL 2h =L 20:5h Þ= log 2. The L 2 -norm curve of each case is also given in Fig. 8.It can be found that the slope of the L 2 -norm curve of the water-only case is over 2.
For the lossy-medium only case, even though the FCD is a second-order method (Celik and Duman, 2012;Liu et al., 2015;Qin et al., 2018), the order of the scheme is below 2.
For the case with the ring-shaped lossy material, the order of the simulation accuracy is back to above 2, but slightly lower than the water-only case.
FIG. 1. (Color online) Grid mesh around an irregular-shaped lossy medium in the Cartesian coordinates surrounded by water.
FIG. 2. (Color online) Pressure contours for plane wave propagation at t ¼ 22:5 ls in (a) the lossy medium, (b) water.

TABLE I .
L 2 -norm pressure errors and observed orders of accuracy for water-only, lossy medium only, and ring-shaped lossy medium cases.
2 -norm at grid size (m) L 2 -norm at grid size (m) L 2 -norm at grid size (m) L 2 -norm at grid size (m) Observed