Accurate simulation of transcranial ultrasound propagation for ultrasonic neuromodulation and stimulation Open Access

The use of implanted electrodes for deep brain stimulation (DBS) is a well-established, invasive treatment for multiple neurological conditions and has directly resulted in a greater understanding of functional neuroanatomy and deep brain circuitry.¹ Unfortunately, its usefulness is limited by the inherent risks of the required neurosurgery combined with difficulties in targeting and repositioning the stimulatory focus.² Non-invasive alternatives such as transcranial magnetic and direct current stimulation have both met with success in research and clinical settings. However, they are limited in terms of their ability to achieve tight spatial focusing, and their penetration deep into tissue.³ Table I demonstrates a selection of existing and planned DBS target structures alongside their approximate dimensions and deviation from the approximate center of the brain—the mid-commissural point (MCP).^4,5 These dimensions demonstrate the millimeter scale size of the target structures, and their position close to the center of the brain. Thus, the ability to non-invasively target these nuclei for modulation and stimulation would represent a revolutionary neuroscientific tool with both clinical and research applications.

Ultrasonic neuromodulation and stimulation (UNMS) offers a potential solution to these requirements, and has recently received a great deal of interest. Transcranial focusing of ultrasound offers the potential for reversible, non-invasive neural excitation and modulation, with focusing on the scale of the acoustic wavelength.³ Table II shows a selection of UNMS papers published in the last decade, and demonstrates the variety of acoustic intensities and frequencies used, target structures sonicated, and neural responses observed. The physical mechanism underlying UNMS remains unclear, although a non-thermal mechanism is suspected, and lower acoustic frequencies have been shown to evoke a response more reliably.^3,6 Most recently ultrasound has been used to elicit electro-encephalogram (EEG) and sensory responses in human subjects, although this has been restricted to superficial cortical brain areas using unfocused single element transducers.^7–9 If UNMS is to develop as an effective non-invasive neurostimulation technique, its application to human subjects must be extended to deep brain targets. Based on the dimensions of DBS targets shown in Table I, and the range of effective ultrasonic intensities shown in Table II, the following focusing requirements may be defined:

• A spatial targeting error of less than 1.5 mm.

• Control of the intensity at the focus with ≤10% error will ensure that neurostimulation occurs. Greater accuracy may be desirable in studies of the mechanisms and thresholds of UNMS.

• An ultrasonic stimulation focus of no greater than 3 mm diameter will ensure stimulatory specificity.

• Steering of the ultrasonic focus up to ∼30 mm from the MCP to allow stimulation of arbitrary deep brain targets.

The primary obstacle to achieving these ultrasonic focusing criteria within the brain is the presence of the skull, which aberrates and attenuates incoming wavefronts. Time-reversal (TR) focusing, first adapted for transcranial focusing by Aubry et al., is able to correct for the aberrating effect of the skull by taking advantage of the time-symmetry of the lossless acoustic wave equation.⁹ In model-driven TR, numerical models simulate the propagation of ultrasound from a target area to a virtual transducer using acoustic property maps of the head derived from CT or MRI images.^9,10 The pressure time series at the simulated transducer surface is then time-reversed, and used to generate drive signals for a multi-element acoustic transducer array. For high-intensity thermal applications, model-driven TR may be combined with MRI thermometry for treatment verification. Chauvet et al.¹¹ confirmed the potential for model-driven TR-based focusing inside the human head to millimeter precision, verified by MRI thermometry. Marquet et al.¹² showed that model-driven TR is capable of restoring 90% of the peak pressure that can be obtained with gold-standard hydrophone based methods when focusing through an ex vivo skull bone. However, model-driven TR remains subject to systematic errors and uncertainties with a resulting loss in focusing quality or efficiency. Four categories of uncertainty are:

• The underlying physical model and how the governing equations model the physics of propagation including phenomena such as absorption, nonlinearity, and shear wave effects.

• Numerical approximations due to the discretization of the simulation domain, including numerical dispersion, the representation of medium heterogeneities, and the effectiveness of any absorbing boundary conditions.

• The inputs to the model, such as the map of acoustic medium properties and the representation of acoustic transducers.

• How the numerical simulations are used within a broader TR protocol, including how the simulated source is related to the desired pressure at the target, and how phenomena that are not time-invariant, such as absorption, are accounted for.

TR simulations for transcranial focusing have typically made use of finite-difference time domain (FDTD) numerical models.^11,12 Recently a k-space corrected, pseudospectral time domain (PSTD) numerical scheme was used in model-driven TR and shown to give comparable accuracy.¹³ Both FDTD and PSTD schemes use consistent approximations to the wave equation, and can be made stable by choosing the discretization parameters appropriately. As the rate of spatial and temporal sampling increases, they will converge on the true solution at a rate dependent on the particular approximations of the numerical model [(ii) above]. However, due to the large scale of these simulations, it is desirable to minimize the grid size and resulting computational burden without compromising accuracy, so knowledge of the minimum sampling criteria necessary to achieve the required accuracy is valuable. In the present paper, these numerical schemes are briefly described, and the various factors affecting the rate of numerical convergence are examined. This is quantified in terms of the spatial and temporal sampling required to obtain acceptable accuracy in the simulation of ultrasound propagation from the scalp to a deep brain target. While the criteria used are established for the application of transcranial UNMS, these results are also applicable to other therapies that require accurate transcranial ultrasound simulation, such as high intensity focused ultrasound (HIFU) ablation and opening the blood brain barrier with ultrasound.

FDTD methods have seen extensive use in the simulation of ultrasound propagation, and have accordingly been used for the purpose of model-driven TR with success.^9–12 In finite difference methods, partial derivatives are calculated using a linear combination of function values at neighboring grid points. The finite difference approximations are derived by combining local Taylor series expansions truncated to a fixed number of terms.¹⁴ When simulating ultrasound propagation, this approximation causes an unphysical dependence of the simulated sound speed on the number of grid points per wavelength (⁠ $PPW = λ / Δ x$⁠) and the number of temporal points per period [ $PPP = 1 / ( f Δ t )$] where $f$ and $λ$ are acoustic frequency and wavelength, respectively, and $Δ x$ and $Δ t$ are the spatial and temporal discretization, respectively.¹⁴ This manifests as a cumulative error in the phase of propagating waves, termed numerical dispersion. In addition, stability conditions must also be met to ensure the numerical scheme is stable. These conditions are contingent on the exact scheme used and the number of simulated dimensions.¹⁴ A useful metric when considering stability is the Courant-Friedreichs-Lewy (CFL) number, defined as

where $c$ is the sound speed. Stability criteria are often expressed as limits placed on the CFL number.^14–16 

In PSTD methods, spatial derivatives are calculated by decomposing the spatially varying acoustic variables into a finite sum of weighted global basis functions.¹⁷ This decomposition allows efficient computation of spatial derivatives using the derivatives of the basis functions. For wave problems, a Fourier basis is typically used, with the basis function weights calculated via the fast Fourier transform.¹⁷ The subsequent gradient calculation is exact, eliminating numerical dispersion due to spatial discretization. However, for an explicit time-stepping scheme, temporal gradients must still be approximated via a finite difference method, with resulting dispersive error.¹⁵ Fortunately, for a second-order accurate approximation, this error can be calculated analytically, and used to introduce a correction factor, $κ = sinc ( c ref k Δ t / 2 ) ,$ in the spatial frequency domain.¹⁵ Here, $k = k x 2 + k y 2 + k z 2$ is the magnitude of the wavevector $( k x , k y , k z )$ at each grid point in the spatial frequency domain, and $c ref$ is a user defined reference sound speed. This method is often called the k-space PSTD method, and in homogeneous media it is unconditionally stable and free from numerical dispersion for arbitrarily large $Δ t .$^15,16 In media with a heterogeneous sound speed, the application of $κ$ in the spatial frequency domain means that a single sound speed must be chosen for the correction factor. As a result, numerical dispersion will arise, the extent of which will depend on the temporal sampling and the difference between the local sound speed $c ( x ) ,$ and the reference sound speed $c ref .$¹⁶ As with FDTD schemes, simulation-dependent limits on the CFL number must be observed to ensure numerical stability.

FDTD and PSTD methods both use functions to interpolate between the values of the acoustic variables at the grid points. The interpolating functions are used to approximate the field gradients at these points, and the values of the field variables are updated at the grid points at each time step. FDTD methods use polynomials to interpolate between neighboring points, while PSTD methods use a Fourier series to interpolate between all points simultaneously.¹⁷ This Fourier series is bandlimited (truncated) to ensure a unique Fourier representation and is therefore known as the bandlimited interpolant (BLI).¹⁷ This can be considered the representation of the discretely sampled pressure field within PSTD schemes. When a time-varying source is used, the resulting pressure signal is formed from a sum of one or more weighted BLIs. As a result of this, a discrepancy can arise between the BLI and the intuitive expectation of what the sampled function represents. This is shown in Fig. 1(a) for a Kronecker delta represented on a discrete grid. In this case, because the Fourier coefficients of the sampled function do not decay to zero before the Nyquist limit of the grid, the intended field is replaced with a BLI representation with Gibbs type oscillations. It is important to understand that this representation is not erroneous per se, but that there is a disparity between the desired input to the PSTD scheme (in this case a Kronecker delta), and what the scheme is capable of representing via a bandlimited Fourier series. To reduce the size of the disparity, smoothing of the intended field can be used to force the Fourier coefficients to decay.¹⁸ This is shown in Fig. 1(b) for the same Kronecker delta function when frequency is filtered with a Blackman window. Although this remains an inexact representation of the original Kronecker delta function, the non-oscillating BLI more closely matches the intended underlying pressure distribution as defined by the values at the discrete grid points.

In this section, the impact of various factors which affect the convergence of FDTD and PSTD models for the case of transcranial ultrasound simulation is presented. These comprise: the influence of the BLI, changes in the effectiveness of the absorbing perfectly matched layer (PML), the impact of numerical dispersion, the representation of discontinuities in medium properties, and staircasing of acoustic sources and material geometry. The first two represent fundamental considerations in numerical simulations, and are dealt with independently. For the subsequent phenomena, the specific inaccuracies occurring when simulating the propagation of ultrasound from a source in the deep brain to an external transducer are established. This is modelled as consisting of 10 cm propagation through cerebral soft tissue, 1 cm propagation through bone, and 1 cm additional propagation through superficial soft tissue, shown in Fig. 2. Accuracy is quantified in terms of the resulting error in the amplitude and position (calculated using time of arrival) of the temporal maximum intensity at the target position and the sampling criteria constraining these errors below 10% and 1.5 mm, respectively, are established. Beam steering capabilities are determined primarily by hardware, and are not considered here. Modeling of the skull as a single homogeneous layer in this way was recently validated for low frequency model-driven TR by Miller et al.¹⁰ The influence of reverberations within the head is considered when necessary, with each reverberation consisting of 2 cm propagation through bone, and 20 cm propagation through brain tissue. The combined effects of these numerical inaccuracies and the validity of the established sampling criteria are then examined through convergence testing of fully simulated two-dimensional (2D) and three-dimensional (3D) TR protocols.

Numerical simulation of ultrasound was carried out with the open source k-Wave toolbox using PSTD and k-space corrected PSTD (with a user defined $c ref$⁠) numerical schemes.^16,19 These are henceforth referred to as “PSTD” and “k-space” schemes, respectively. The toolbox also includes a second-order accurate in time, fourth-order accurate in space (2–4) FDTD scheme, which was also tested. This scheme is described in detail by Strikwerda,¹⁴ and is widely used to simulate acoustic wave propagation, including in simulated TR.^11,12 Unless stated otherwise, the CFL number was set to 0.3, one-dimensional (1D) tests were carried out on a spatial grid of 4096 grid points, and 2D tests on a 1024 $×$ 1024 grid. Frequency filtered Kronecker delta functions, like that shown in Fig. 1(b), were used to create broadband pressure sources. Homogeneous simulation grids were given the acoustic properties of brain tissue, a density of 1040 kg/m³, and a sound speed of 1560 m/s (also used to represent superficial soft tissues).²⁰ For heterogeneous simulations, bone tissue was assigned a density of 1990 kg/m³ and a sound speed of 3200 m/s.²⁰ When it was necessary to define a specific ultrasonic frequency of interest to calculate the required sampling criteria, 500 kHz was used. This frequency has seen extensive use in studies of UNMS (see Table II), sits within the range of ultrasound frequencies demonstrating optimal transcranial transmission,^21,22 and has a theoretical minimum focus size of ∼3 mm diameter in soft tissue.

The BLI represents a fundamental component of both k-space and PSTD schemes. As such, it is necessary to examine its impact on simulation accuracy before moving on to more complex factors that affect the rate of convergence. Bandlimited interpolation, as described above, can result in a discrepancy between the intended pressure field and the representation of that field within PSTD schemes when the Fourier coefficients of the intended field have not decayed sufficiently. Practically, this manifests globally as undesired, oscillating pressure values across the simulation grid [see Fig. 1(a)]. Therefore, to examine the impact of BLI effects, it is necessary to determine the amplitude of these undesired pressures relative to that of an intended input.

In practice, the error in the representation of a particular pressure distribution will depend on how well it can be represented by a discrete Fourier transform at a specific spatial discretization.¹⁷ Tonebursts have a well-defined power spectrum determined by their length and central frequency. Therefore, to approximate the BLI error likely to be generally observed, a series of time-varying 10, 30, and 50 cycle acoustic toneburst sources with central wavenumbers approaching the spatial Nyquist limit were used as input signals. These sources have 22.7%, 7.4%, and 4.3% full width at half maximum (FWHM) bandwidth as a percentage of central frequency, respectively. The source was positioned a quarter of the way along a homogeneous 1D computational grid with no PML. The simulations were run for the time taken for waves to travel from the source to the center of the grid. The pressure was recorded at every grid point of the other half of the grid which, according to causality, should have remained quiescent if the BLI of the pressure field matched the intended input of compactly supported tonebursts. Error was quantified as the maximum pressure recorded across the second half of the grid relative to the peak pressure of the source toneburst. The results are shown in Fig. 3(a) as a function of the PPW of the central wavenumber of the toneburst. The amplitude of the non-causal pressure drops rapidly as the number of PPW increases from the Nyquist limit. Reducing the error requires a higher number of PPW for shorter tonebursts due to their wider power spectra, but for all three toneburst lengths the error drops to below −60 dB by 3 PPW. An additional observation was that wavenumbers corresponding to less than 2 PPW are not aliased or otherwise propagated on the grid.

To determine what frequencies comprised the observed non-causal pressure, the results obtained from 10 cycle tonebursts were examined further. Time-varying pressure signals were recovered from the grid points closest to the wave front, which experienced the peak non-causal pressures. The normalized amplitude spectra of these signals resulting from source tonebursts with central wavenumbers sampled at 2 and 2.4 PPW are displayed in Fig. 3(b), alongside the corresponding amplitude spectra of the source tonebursts. The recorded spectra demonstrate a sharp peak at 2 PPW regardless of the central frequency of the source toneburst, and the amplitude of the peak scales with the amplitude of the 2 PPW component of the source toneburst. Practically, these results demonstrate that this error reduces very rapidly as the spatial sampling of the pressure distribution increases, and at 3 PPW BLI errors are reduced to below −60 dB. In higher dimensions, BLI errors are less severe than the 1D case.

The interaction of outgoing pressure waves with the edge of the simulation grid presents a problem for numerical schemes. In the FDTD scheme used here, outgoing pressure waves are perfectly reflected from the edge of the grid, while for k-space and PSTD schemes outgoing pressure waves are “wrapped” to the opposite edge of the simulation grid (i.e., the grid is toroidal). To replicate free field conditions, k-Wave employs Berenger's split field PML, where the pressure field is artificially divided into Cartesian components to allow selective absorption of the normally incident component.^19,23 The response of the PML to different frequencies was established by propagating broadband pressure sources toward a 20 point PML on a homogeneous 1D grid using the PML profile given in Tabei et al.^15,19 Rather than the 2-4 FDTD scheme used elsewhere, a second order accurate in space and time (2-2) FDTD scheme with a CFL of 1 was used to prevent numerical dispersion. It should be noted that this FDTD formulation is not practical outside the homogeneous 1D case due to stability constraints. The time-varying pressure traces of the incident wave, the wave reflected from the surface of the PML, and the wave transmitted to the edge of the computational grid were recorded. The power spectra of these signals were used to calculate reflection and transmission amplitudes relative to the incident wave as a function of spatial sampling. No notable difference was observed between k-space and PSTD schemes. Results are shown in Fig. 4 for k-space and 2-2 FDTD schemes. In both schemes, the pressure reflection coefficient demonstrates a dependence on spatial sampling, rising steadily from below −120 dB for frequencies sampled at above 4 PPW to total reflection at 2 PPW. Transmission to the edge of the grid remains constant at below −70 dB for both schemes until spatial sampling drops beneath 3 PPW, below which the k-space scheme shows an increase in transmission and the FDTD scheme shows a reduction in transmission. These results indicate that the effectiveness of the PML is greatly reduced for wavenumbers sampled at below 3 PPW, and it cannot be relied on at these PPW values. However, erroneous reflection and transmission reduce rapidly as sampling increases. It should be noted that pressure reaching the edge of the grid for both schemes is subject to further attenuation within the PML when reflected or wrapped back into the grid. Furthermore, the BLI will have influenced the behavior of these tests for frequencies sampled at close to the Nyquist limit, which may explain why the k-space scheme shows an increase in both reflection and transmission close to 2 PPW.

The PML was also tested in 2D to determine its dependence on the angle of incidence of incoming waves. A broadband point source was placed close to the edge of the PML on the 2D grid and propagated into, and across the surface of, the PML. The pressure was recorded at the edge of the simulated domain to examine transmission, with each recording position corresponding to a particular angle of incidence. The peak pressure transmission at each angle was calculated through comparison with a reference recording obtained with PML absorption set to zero. The results are shown in Fig. 4(c). Transmission to the edge of the grid is lowest for normally incident waves, rising with increasing angle of incidence crossing to above −60 dB at ∼40°. No clear relationship between angle of incidence and reflection from the PML was observed. These results should be considered when designing acoustic sources and considering the angles at which pressure waves will impinge on the PML.

To examine the impact of numerical dispersion, a broadband pressure source was defined on a homogeneous 1D grid. Grids with the medium properties of both bone and brain tissue grids were tested, and for the k-space scheme c_ref was set to the speed of sound in brain tissue. For FDTD schemes, the temporal and spatial dispersive errors oppose each other, with reduced dispersive error at higher CFL numbers.¹⁴ Therefore the CFL was set to 0.5 for these simulations, the highest value at which both schemes are stable in 3D.^14,15 The time-varying pressure was recorded at a distance of 1 cm, and the phase spectra of the recorded pressure signals were compared to a dispersion-free reference simulation obtained with perfect k-space correction. This allowed calculation of phase error per cm propagated in either tissue type as a function of acoustic frequency. Using the model for transcranial propagation of ultrasound to a deep brain target described in Sec. III A, this was used to calculate the sampling criteria required to obtain <1.5 mm positional error for the direct path, and for each reverberation, shown in Table III. To compare, when sampled at 18 PPP, the k-space scheme is exact for soft tissue and gives a 19 μm error in the focal position per cm propagated in bone, PSTD gives a 25 μm error per cm in soft-tissue and a 50 μm error per cm in bone, and FDTD gives 25 μm error per cm in soft-tissue and a 130 μm error per cm in bone.

Results are given in terms of temporal PPP to allow a comparison between dispersive error for the FDTD scheme, which is dependent on both spatial and temporal sampling, and PSTD and k-space schemes, which are dependent on temporal sampling only. Equation (1) shows how the CFL defines a fixed ratio between spatial and temporal sampling. Different CFL numbers will result in a different combination of spatial and temporal requirements for the FDTD scheme. Values for both PSTD schemes are dependent only on temporal sampling, and do not require a particular spatial sampling to reduce dispersive error. However, with that in mind, the results shown in Table III do demonstrate a clear reduction in the temporal sampling required to minimize dispersive positional error below acceptable levels for the k-space scheme compared to FDTD and PSTD schemes.

To examine the delay in convergence due to inaccuracies in reflection and transmission from medium discontinuities, broadband pressure sources were propagated across a bone-soft tissue interface (propagation direction makes no difference). The incident, reflected, and transmitted waves were recorded and the power spectra used to calculate intensity reflection and transmission coefficients for each wavenumber. Percentage error in these coefficients was calculated through comparison with the analytical values. To examine the dependence of this error on the size of the impedance change, these tests were repeated with the impedance of the bone varied up to ten times that of the soft tissue, with sound speed and density varied independently. No difference was observed between k-space and PSTD schemes.

Figure 5 shows the resulting error in intensity transmission and reflection coefficients as a function of PPW, and as a function of impedance change for a PPW of 6. FDTD and k-space schemes demonstrate a similar error even at high impedance changes. Changes in density result in increased deviation from the correct coefficient due to the interpolation of the density values on a staggered grid.

To calculate how the error in 1D reflection and transmission will affect transcranial focusing, the error in the intensity following transmission across a bone layer (i.e., across two bone/soft-tissue interfaces) was computed as

where $Te$ is the analytical energy transmission coefficient between bone and soft tissue, and $Te ̃$ is the simulated energy transmission coefficient as a function of spatial PPW. Reflections inside the skull and skull cavity were not considered. To obtain <10% intensity error, the FDTD scheme requires 5.9 PPW while the k-space and PSTD schemes require 4.3 PPW. This result is notable, as the representation of discontinuities in medium properties has previously been identified as a limitation of PSTD methods.¹⁵ During the update steps of these schemes, the pressure field is multiplied by the maps of medium density and sound speed, before being evaluated by a truncated Fourier series. Step changes in medium properties will therefore introduce Gibbs phenomenon into the pressure field, as described in Sec. II C. However, these results indicate that, for the step change in medium properties between bone and soft tissue, the error resulting from the representation of this change within the FDTD scheme tested is greater.

Staircasing refers to the spatial approximation that is necessary when attempting to define continuous geometries on a discrete, regular Cartesian grid in 2D and 3D. Curved surfaces and lines at an angle to the Cartesian directions will be approximated in a stair-stepped manner, and certain vertex and edge positions do not correspond to points on the grid.²⁴ The impact of staircasing was examined separately for acoustic sources and medium distributions. Tests involved recording the time-varying pressure at a number of positions across the field resulting from a staircased representation of a source or medium, and comparison of these signals with references obtained from a staircase free simulation. Error was then quantified as the percentage error in the amplitude of the temporal peak intensity (calculated using a plane wave assumption) and its positional error (derived from the change in the time of arrival of the intensity peak) as a percentage of wavelength. A positional error of 50% of wavelength corresponds to 1.5 mm for a source frequency of 500 kHz in brain tissue. These errors were calculated for each recording position, and then averaged across the field to give mean errors in peak intensity amplitude and position. No notable difference in error was observed between FDTD, PSTD, and k-space schemes across all tests.

The impact of staircasing on acoustic sources was examined using line-sources with a length of $65 × d x$⁠, where $d x$ is the spatial discretization step, at a series of angles to the Cartesian grid. These included four angles that form Pythagorean triangles on the grid, specifically: ∼14.3° (with Pythagorean triple 16, 63, 65), ∼22.6° (25, 60, 65), ∼30.5° (33, 56, 65), and ∼36.9° (39, 52, 65). For these angles, the line-source endpoints are coincident with specific grid point positions, and any error is only due to the staircased representation of the line, rather than endpoint misregistration (both are aspects of staircasing error). A source defined parallel to a Cartesian axis was used as a non-staircased reference. The sources were excited with 10 cycle acoustic tonebursts with a range of central wavenumbers sampled at 3–100 PPW. The amplitudes of the source signals were normalized based on any change in the number of distinct source points used when defining an angled line source when compared to the aligned case. The time-varying field was recorded at 100 points positioned in front of the line source, and the sensor map was rotated with the line source to maintain source-sensor geometry. The simulation layout is shown in Figs. 6(a) and 6(b).

Mean errors in the amplitude and position of the peak intensity across the sensor field were calculated independently for each angle tested relative to the aligned, non-staircased reference case. The maximum mean errors across the range of angles tested for each PPW value are shown in Fig. 7. These results demonstrate that staircasing errors worsen with lower spatial sampling and are less serious for Pythagorean angles, when endpoints are correctly registered. The error in the position of the intensity peak never rises above 50% of wavelength for any source. Seventeen PPW are required to obtain <10% error in the amplitude of the intensity peak for all angles tested, while Pythagorean angles require only 7 PPW. Although the error examined here does not relate directly to the model for transcranial ultrasound propagation described above, these results do indicate that staircasing and spatial sampling must be considered when defining acoustic source distributions, and that error can be reduced by ensuring endpoint registration. The testing of multiple angles also allowed examination of how the exact mapping of the staircased line relates to the error observed in the resulting field. No clear relationship between the angle of the line source and the level of staircasing error was observed. However, a staircasing metric was defined as the average distance between the staircased source points and their equivalent equispaced points on an ideal angled line. It was observed that the convergence rate of the error of the staircased source maps, quantified as the average L2 error in the recorded pressure signals at the maximum spatial sampling tested, showed a strong dependence on this staircasing metric. Although only a simple metric, this indicates that the severity of staircasing error can be predicted through comparison of an ideal or parametric map of the intended geometry with its staircased representation.

The impact of staircasing of heterogeneous medium properties was examined in two separate tests. The first was conceptually similar to the examination of source staircasing. An acoustic point source excited by ten cycle acoustic tonebursts with central frequencies ranging from 3 to 80 PPW was propagated across a planar medium boundary (soft tissue-to-bone), defined at varying angles to the Cartesian axes. The time varying pressure field was recorded at 100 sensor points following the interaction with the medium boundary. The sensor points were rotated with the medium boundary to maintain the simulation geometry. A non-staircased boundary defined along a Cartesian axis was used as a reference map. This simulation layout is shown in Figs. 6(c) and 6(d). The second test was designed as a more accurate model of staircasing in transcranial transmission. A 10 cycle, 2520 PPW toneburst was propagated through a medium map comprising a quarter circle bone layer. The medium was then artificially staircased through spatial downsampling, before being remapped to the original grid. A 3780 $×$ 3780 simulation grid was used due to the large number of integer factors of 3780, which allowed the medium to be successively downsampled while maintaining positioning. The time-varying pressure was recorded across a quarter-circle, and error metrics computed through comparison with the least staircased medium distribution. This simulation layout is shown in Figs. 6(e) and 6(f). An effective PPW value for each level of downsampling was calculated through comparison of the source PPW with the new effective spatial discretization. Mean error measurements across the recorded fields as a function of PPW are shown in Fig. 8 for both medium staircasing tests. For the single interface model, the maximum mean errors across the range of angles tested are shown.

In terms of the impact on the model for transcranial propagation, the results for the bone layer model shown in Fig. 8 suggest that 20 PPW or above are required to obtain less than 10% mean error in intensity transmitted to an external transducer surface. As might be expected, the errors for the single interface are lower, with 6 PPW required to obtain the same error in peak intensity amplitude. The error in peak intensity position is less serious, with mean positional error never rising above 50% of wavelength (1.5 mm for 500 kHz ultrasound in brain tissue) for both tests, as with source staircasing. This may be due to staircasing introducing a random error in acoustic pathlength, leading to a defocusing and change in amplitude rather than a shifting in the peak position. To place this in context, the voxel size of clinical CT images is on the order of 0.5 mm at best. This corresponds to 6 PPW at 500 kHz, and fewer at higher frequencies. This suggests that staircasing may have a significant impact on simulations using image derived medium property maps. When considered alongside the results for source staircasing, these results indicate that staircasing error is likely the most serious of the numerical errors tested. The single interface medium model was also briefly tested using an elastic PSTD model, which indicated that medium staircasing may also have a pronounced impact on simulated mode conversion, although more rigorous testing is necessary.²⁵ 

To examine the combined effects of numerical errors on the effectiveness of transcranial TR focusing, convergence testing of fully simulated TR was carried out in 2D and 3D. The general method is outlined in Fig. 9. It consists of forward propagation of a 10 cycle toneburst from a source point inside a virtual skull model to a circular (2D) or hemispherical (3D) virtual transducer-sensor array. As these tests were to examine the specific impact of numerical convergence, the simulated transducer array surface was modelled as a continuous surface made up of point transducers at the resolution of the spatial grid, and no attempt was made to replicate real transducer characteristics. This ensures that convergence is dependent only on the numerical accuracy of the forward simulations and will apply to alternate source conditions. The spatial discretization of these simulations was varied to correspond to a range of spatial PPW values for the central frequency of the source toneburst. The time-varying pressure signals recorded at the virtual transducer position were then time reversed and propagated into the head to refocus onto the target position. The reversal simulations were carried out at the finest discretization feasible, and the CFL was 0.3 for all simulations. Due to the change in spatial and temporal discretization, it was necessary to interpolate the pressure signals recorded in the forward simulations onto the spatial and temporal grids used in the reversal simulations [see Fig. 9(b)] using Cartesian triangulation and Fourier interpolation, respectively. In addition, the position of the source, defined at a grid point on the high resolution reversal grid, was not conserved due to the varying discretization of the forward simulations, and instead the nearest neighboring point was used. The peak pressure occurring in a time window of 20 acoustic cycles centered on the expected time of refocusing was recorded across the brain volume. Convergence was established by examining focusing quality as a function of the discretization used in the forward simulations. The focusing metrics examined were the spatial and temporal peak pressure across the brain volume, the distance of the peak from the target position, and the FWHM of the focal spot size. Peak pressure and focus FWHM were normalized relative to the results obtained for the most highly resolved forward simulation.

Simulations in 3D were carried out with toneburst sources with central frequencies of 500 kHz, while testing in 2D used both 250 and 500 kHz tonebursts. The simulation parameters employed across these tests are shown in Table IV. Grid sizes include the absorbing PML layer. The forward simulations were run for the time taken for an acoustic wave to propagate across the grid three times, plus the duration of the source toneburst. The reversal simulation was run for the additional time of five acoustic cycles in order to fully capture the reconstructed toneburst. The medium property map used in the convergence tests was derived from a T1-weighted MR image obtained from the Imperial College brain development dataset.²⁶ Brain and skull volumes were extracted using the FSL MRI processing toolbox²⁷ and converted into a surface mesh using the iso2mesh toolbox.²⁸ The reference surface mesh was then sampled onto a 2D or 3D grid with the required spatial discretization steps. Examples of the 2D maps used in forward and reversal simulations, and their varying discretizations, are shown in Figs. 9(a) and 9(c). In each case, the skull was modelled as a single homogeneous bone layer, and the rest of the simulation domain was assigned brain tissue medium properties. Although fully heterogeneous models of the skull have demonstrated tighter model-driven TR focusing in some cases, homogeneous models remain effective.²⁹ Furthermore, they allow accurate resampling of the bone map to multiple spatial discretizations without interpolation, and ensure that convergence is dependent on the accuracy of the numerical simulation, rather than the mapping and homogenization of acoustic medium properties. For the k-space scheme $c ref$ was set to the speed of sound in brain tissue.

2D convergence testing was carried out using 10 cycle acoustic toneburst sources corresponding to 250 and 500 kHz frequencies. Reversal simulations were carried out using a 3072  $×$ 3072 grid, with a spatial discretization corresponding to 50 PPW for 500 kHz and 101 PPW for 250 kHz. Forward simulations were carried out using the k-space and FDTD schemes. Reversal simulations were carried out using the k-space scheme only.

The results of the 2D convergence testing are shown in Fig. 10, with error in peak pressure position given relative to the source point in the forward simulations. Refocusing quality in the reversal simulations increases with the spatial PPW of the simulated frequency in the forward simulations. Several key points can be derived from these results. First, for all three metrics, the k-space scheme demonstrates convergence at approximately ∼2 PPW below the FDTD scheme. Second, although there is some difference in the position and size of the focus at very low sampling, both 250 and 500 kHz demonstrate similar behavior as a function of spatial sampling. This indicates that these results can, to some extent, be generalized, and suggests that the reversal simulations have converged for both frequencies. Finally, of the three refocusing metrics examined, normalized peak pressure across the brain volume [Fig. 10(a)] requires higher spatial sampling to converge than either focal volume or the deviation of the focus from the target. This suggests that when fine pressure control is not required, coarser sampling criteria may suffice.

3D convergence testing employed a similar protocol to the 2D convergence testing described above. Ten cycle acoustic tonebursts with a central frequency of 500 kHz were propagated from a source point inside a full 3D model of the human skull to a simulated hemispheric transducer. Reversal simulations were carried out using a 1024  $×$ 1024 $×$ 1024 simulation grid with a spatial discretization corresponding to 16.7 PPW. Forward and reversal simulations were carried out using the k-space scheme only. Testing was also carried out using homogeneous media for forward and reverse simulations, to test the accuracy of the spatiotemporal interpolation. The results for heterogeneous 3D convergence testing are shown in Fig. 11. The results for normalized peak pressure amplitude in Fig. 11(a) show similar trends to the 2D results. Six PPW are required to obtain 95% reconstruction of pressure at the target (corresponding to <10% drop in intensity) and ∼10 PPW to attain convergence. The convergence of the volume of the focal spot [Fig. 11(c)] shows similar behavior, although it only requires ∼6 PPW to fully converge. Given the known impact of staircasing in 2D, the faster convergence here is likely due to a reduced staircasing error for 3D geometry. The results in Fig. 11(b) show the error in the position of the pressure peak relative to both the reversal target and the shifted forward source point. This error is reduced compared to the 2D case, never rising above 1.1 mm and, when computed relative to the position of the forward source, is stable using sampling as low as 2 PPW. However it can be seen from the other results that at this sampling the peak pressure amplitude is much lower, with a larger focal spot. The apparent periodicity in the positional error when calculated relative to the parametrically defined target point is likely due to the oscillating distance of a definable source point from this position. Accordingly, the reduction in error and lack of dependence on spatial sampling when calculating positional error relative to the actual source position from the forward simulation suggests that any error in the position of the peak is in fact due solely to the misregistration of the source points. Generally, these results confirm that when targeting accuracy is the prime concern, spatial sampling requirements are laxer than when a tight focal volume with known peak pressure amplitude is required. This is in agreement with previous studies which have demonstrated that good spatial targeting of HIFU can be obtained via simulated TR using relatively coarse spatial sampling.^11,30 Homogeneous testing demonstrated total convergence across all metrics by 3.5 PPW, which is to be expected given the behavior of the PML and BLI, discussed in Sec. III.

In this paper, a comprehensive assessment of the impact of different factors that affect the convergence of numerical models for the simulation of transcranial ultrasound propagation was carried out. The spatial and/or temporal sampling required to reduce inaccuracies below the levels required for targeting of deep brain nuclei for neurostimulation were determined.

Initial simulations examined reduction in the effectiveness of the PML, and the impact of the BLI when using k-space and PSTD methods. Both the PML and BLI lead to erroneous pressures appearing on the grid when simulating frequencies sampled at close to the spatial Nyquist limit. Although both of these effects have the potential to seriously reduce the accuracy of the simulations, they decrease in severity rapidly as the rate of spatial sampling increases. Above ∼3 PPW, erroneous pressures resulting from both BLI and PML effects were at least −60 dB below the amplitudes of the ultrasound sources being simulated.

Numerical dispersion has a serious effect on the accuracy of FDTD and PSTD schemes, resulting in high temporal sampling requirements to reduce positional error. However, this was not the case for the k-space scheme, where ∼3 PPW will serve to limit dispersion sufficiently for transcranial transmission for any stable CFL value. Errors in reflection and transmission from discontinuous medium properties manifest in the magnitude of reflected and transmitted simulated intensities. Despite the representation of step changes in media previously being identified as a key limitation of PSTD schemes, the error was shown to be more severe for the 2-4 FDTD scheme tested. To reduce error in the intensity below 10% following transcranial transmission, k-space and PSTD schemes require 4.3 PPW, while FDTD requires 5.9 PPW.

Staircasing of source and medium geometries was shown to require the most stringent sampling criteria to obtain required accuracy, affecting FDTD, PSTD, and k-space schemes equally. Both source and medium staircasing were shown to have a greater impact on the intensity amplitude of the toneburst signal being examined than the position of the intensity peak. The results shown in Figs. 7 and 8 indicate that $≥$20 PPW are required to reduce the error in peak intensity following transcranial transmission below 10%. The preliminary examination of a potential staircasing metric also suggests that the error resulting from a particular staircased geometry is directly related to its deviation from the ideal geometry.

Convergence testing of a fully simulated TR protocol using 2D and 3D head models was used to examine the impact of all numerical errors in concert. Testing in 2D for 250 and 500 kHz ultrasound showed a faster rate of convergence for all focusing metrics for the k-space scheme when compared to FDTD. In addition, the error in the peak pressure amplitude at the focus showed slower convergence than both the positional error, and the volume of the focus. This is likely due to the most serious source of error, medium staircasing, which was shown to have a greater impact on the peak intensity amplitude of transcranially transmitted ultrasound, than the position of the peak. Results in 3D showed similar trends to the 2D results for the convergence of the peak pressure amplitude. The focal spot size showed slightly slower convergence in the 3D case, while the positional error demonstrated almost no dependence on the sampling rate of the forward simulation. This indicates that less stringent sampling may suffice for applications concerned only with the position of the focus, rather than the size of the focal spot and the exact amplitude at the target. When fine control over the pressure amplitude is required, stricter sampling may be necessary. Despite the relatively severe error resulting from staircasing at higher spatial sampling, all three metrics of focusing quality were well converged at below 20 PPW. This discrepancy may be due to the differences between the convergence testing protocol and the specific test used to examine staircasing across a bone layer, and suggests that the influence of staircasing is case specific.

The work described above is subject to some limitations, primarily the degree to which the examination of individual numerical errors can be generalized to different setups, although trends and qualitative observations remain valid. Many of the tests only examine toneburst sources, and the error is evaluated over a small field, with pressure recorded at a limited number of sensor positions (see Fig. 6). A separate 2-2 FDTD scheme was used to examine the effectiveness of the PML, which may not be exactly relatable to the commonly used 2-4 FDTD scheme. Furthermore, the impact of shear wave propagation was not examined. This will not have affected 1D or homogeneous simulations, but a more thorough examination of medium staircasing should include testing of elastic wave propagation. Similarly, no effort was made to examine the manifestation of numerical errors when modeling nonlinear propagation or acoustic absorption, which will become relevant for applications requiring the simulation of high-amplitude ultrasound, such as HIFU. It should be noted that, in simulated TR, accounting for acoustic absorption occurs in the post-processing stage, when converting recorded pressure signals into driving amplitudes,¹² and work examining absorption should focus on this stage of the simulated TR process.

The results presented here are primarily relevant to the simulation of transcranial ultrasound propagation for TR targeting of deep brain structures with finely controlled ultrasound for the purposes of neurostimulation. However, the criteria and simulations presented are also relevant to alternative low-intensity, transcranial ultrasonic therapies such as opening the blood-brain barrier with ultrasound,²² as well as existing transcranial HIFU ablation therapies. Use of appropriately discretized simulations will ensure accurate targeting and effective therapy as the field of ultrasonic neurostimulation develops.

The authors would like to thank Charlotte Stagg and Adam Thomas for helpful discussions, and Nishant Ravikumar and Zeike Taylor for the provision of the skull and brain meshes. This work was supported by the Engineering and Physical Sciences Research Council (ESPRC), UK. J.J. is financed from the SoMoPro II programme. This research has acquired a financial grant from the People Programme (Marie Curie Action) of the Seventh Framework Programme of EU according to the REA Grant Agreement No. 291782. The research is further co-financed by the South-Moravian Region. This work reflects only the author's view and the European Union is not liable for any use that may be made of the information contained therein. Computational resources were provided by the IT4Innovations Centre of Excellence project (CZ.1.05/1.1.00/02.0070), funded by the European Regional Development Fund and the national budget of the Czech Republic via the Research and Development for Innovations Operational Programme, as well as Czech Ministry of Education, Youth and Sports via the project Large Research, Development and Innovations Infrastructures (LM2011033).

Simulations were carried out using the open source k-Wave toolbox for Matlab, c++. The toolbox includes k-space, PSTD, and 2-4 FDTD codes for the time-domain simulation of acoustic fields. 1D simulations were carried out in the Matlab environment on a Dell Precision T1700 with an Intel Xeon E3-1240 3.40 GHz CPU and 16 GB of RAM running Windows 10 64 bit. 2D simulations were carried out in the Matlab environment with CUDA hardware acceleration on a Dell PowerEdge R730 compute server with 2 × 6-core Xeon E5-2620 2.4 GHz CPUs, 64 GB of 1866 MHz memory, on an Nvidia Titan X GPU with 3072 CUDA cores and 12 GB of memory. The largest 2D simulations had a domain size of 3780² including the PML and comprised 258 462 time steps, with a total runtime of 10.6 h. 3D simulations were carried out on the IT4I Salomon supercomputing cluster. Each simulation was carried out on Intel Xeon E5-4627v2, 3.3 GHz, 8cores and 256 GB of RAM per simulation. The largest 3D simulations had a domain size of 1024³ including the PML and comprised 22 718 time steps, with a total runtime of 112.3 h. The skull mesh used in convergence testing is Copyright Imperial College of Science, Technology and Medicine 2007. All rights reserved. www.brain-development.org.