Transcranial ultrasound simulations are increasingly used to predict *in situ* exposure parameters for ultrasound therapies in the brain. However, there can be considerable uncertainty in estimating the acoustic medium properties of the skull and brain from computed tomography (CT) images. This paper shows how the resulting uncertainty in the simulated acoustic field can be predicted in a computationally efficient way using linear uncertainty propagation. Results for a representative transcranial simulation using a focused bowl transducer at 500 kHz show good agreement with unbiased uncertainty estimates obtained using Monte Carlo.

## 1. Introduction

Transcranial ultrasound simulations are increasingly used to predict *in situ* exposure parameters for ultrasound therapies in the brain, including for blood-brain barrier opening (Deffieux and Konofagou, 2010), high-intensity focused ultrasound ablation (McDannold , 2019), and transcranial ultrasound stimulation (TUS) (Lee , 2016). Simulations are particularly important for treatment planning in TUS, as online monitoring methods (such as thermometry and acoustic radiation force imaging) are not yet sufficiently sensitive for use with low-intensity ultrasound therapies (Li , 2022).

The current approach for transcranial ultrasound simulation is to map the geometry and acoustic properties of the skull and brain from computed-tomography (CT) images (Marquet , 2009). However, there can be considerable uncertainty in the mapping from Hounsfield units to acoustic parameters (Webb , 2018; Webb , 2021). Errors in the sound speed of the skull, in particular, can result in changes to the amplitude and shape of the acoustic field inside the brain (Montanaro , 2021; Robertson , 2017; Vaughan and Hynynen, 2002). In principle, these variations could be quantified using a Monte Carlo approach, where many simulations are performed using different acoustic property mappings, and the results of these simulations are used to quantify the uncertainty. However, Monte Carlo estimation can take a large number of samples (i.e., simulations) to converge, significantly increasing the computational cost of model-based treatment planning.

In general, quantifying uncertainty is a challenging task for which many methods have been developed to avoid expensive Monte Carlo simulations. These range from the application of the Koopman operator (Gerlach , 2020), to density estimation (Papamakarios and Murray, 2016). Here, we demonstrate how computationally efficient uncertainty estimation can be achieved using linear uncertainty propagation with a differentiable wave simulator, j-Wave (Stanziola , 2022).

## 2. Linear uncertainty propagation

For a given scanner and image acquisition and reconstruction settings, the mapping from a CT image in Hounsfield units to an image of mass density in kg·m^{–3} can be performed with the aid of a calibration image of a test object with known density. Figure 1(a) shows an example of such a calibration acquired using a CIRS electron density phantom (model 062M, Sun Nuclear, Melbourne, FL) using the scanner details described in Caballero-Insaurriaga (2019). The error bars show the standard deviation in the image values averaged over each tissue equivalent electron density plug. A linear curve is then typically used to map from density to sound speed [or, equivalently, two different linear mappings are used to map to density and sound speed from Hounsfield units (Marquet , 2009)]. Figure 1(b) shows a scatterplot of the density and sound speed of the skull samples measured in Webb (2018) and Webb (2021), along with a linear fit and the 95% confidence interval for the linear model. There is considerable spread in the measurement data, which may be attributed to both measurement noise (e.g., in the time-of-flight picking used to estimate the sound speed) and inherent variations in the mapping for different skulls and different regions of the skull. Here, we only consider the uncertainty in the mapping from density to sound speed, but in principle, it is also possible to propagate the uncertainty in the mapping from Hounsfield units to density and to account for uncertainty in other acoustic and thermal parameters.

*ρ*and sound speed

*c*is given by the linear relationship

^{–1}and $ \rho \u0302 = 1782$ kg·m

^{–3}are the mean values for the sound speed and the density in the data, respectively [shown by the cross-hairs in Fig. 1(b)]. Fitting the model on mean-zero data results in a linear fit with uncorrelated estimates for the

*α*and

*β*parameters. While correlation of the input variables can be accounted for by linear uncertainty propagation, this simplifies the implementation. The least squares estimate of the slope and intercept are given by $ \alpha 0 = 1.333$ and $ \beta 0 = 0$ m·s

^{–1}, with a standard deviation of $ \sigma \alpha = 0.168$ and $ \sigma \beta = 18.8$ m·s

^{–1}.

*α*

_{0}and

*β*

_{0}), to account for the uncertainty in the mapping, it is assumed that

*α*and

*β*can be represented by uncorrelated Gaussian probability distributions. For a fixed problem (for example, a particular transducer and density map), the simulated acoustic pressure field

*p*can be written as a function of the parameter vector $ v = ( \alpha , \beta )$,

*f*includes the wave simulation and any other pre- or post-processing steps needed to generate the output pressure field. To propagate uncertainty, it is assumed that $v$ is drawn from a normal distribution centered around

*μ*with variance ∑

_{v}_{v}, given by

*ε*is small enough, Eq. (2) can be expanded as a Taylor series up to second order,

*J*is the Jacobian.

*α*

_{0}and

*β*

_{0}), while the variance is given by Giordano (2016),

*f*at $ v = \mu v$.

## 3. Numerical experiments

To demonstrate uncertainty estimation using Eq. (6), an ultrasound simulation through a CT-derived acoustic property map was performed in the context of TUS. A CT image of the head was taken from the dataset described in Caballero-Insaurriaga (2019) and Miscouridou (2022), re-sampled to an isotropic resolution of 0.5 × 0.5 × 0.5 mm^{3}, and truncated to a 120 × 70 × 70 mm^{3} domain in the region of the primary motor cortex. The simulation uses speed of sound and density values derived from the CT image, while medium is considered lossless. The transducer was a focused bowl defined following Aubry (2022), with a radius of curvature and aperture diameter of 64 mm, source pressure of 60 kPa, and driving frequency of 500 kHz.

Simulations were performed using j-Wave (Stanziola , 2022), a differentiable wave solver, by running a time-domain simulation to steady state and then recording the maximum pressure amplitude over several cycles. To avoid staircasing artefacts, the transducer was represented on the grid using the band limited interpolant (Wise , 2019). The grid spacing gave 6 points per wavelength (PPW) in water, and the time step was set using a Courant–Friedrichs–Lewy number of 0.3.

To calculate the gradients $ \u2202 f / \u2202 \alpha $ and $ \u2202 f / \u2202 \beta $, forward-mode differentiation was used (Stanziola , 2021, 2022). In j-Wave, this can be performed in parallel with the unperturbed forward computation, at a cost of one forward computation for each partial derivative. Forward-mode differentiation, in contrast to the widely used adjoint method, does not necessitate storing the intermediate computations. In the context of wave physics, it has been used, for instance, to perform sensitivity analysis of photonic devices (Hughes , 2019). Alternatively, the gradients could also be approximated using a non-differentiable wave solver by replacing the gradients with a finite difference approximation.

As a reference, the standard deviation of the simulation was also estimated via Monte Carlo, by running 200 simulations (also with j-Wave) with samples of the parameter vector $v$ (i.e., different values for the slope and intercept in the linear relationship between density and sound speed) drawn from its probability distribution, that is, a Gaussian with mean and standard deviation equal to the ones estimated by the linear regression.

## 4. Results

The simulated acoustic pressure field for a plane passing through the focus is shown in Fig. 2(a). The standard deviation calculated using Monte Carlo is displayed in Fig. 2(b), while the standard deviation calculated using linear uncertainty propagation is displayed in Fig. 2(c). The two uncertainty maps are in close agreement, both in terms of the structures shown and the size of the standard deviation estimate. The difference between the two, illustrated in Fig. 2(d), shows that the uncertainty is accurately estimated, particularly below the skull and close to the focal region, while the most noticeable discrepancies between the two are found inside the skull and in the reflected wave-field. The uncertainty estimates for lateral and axial profiles through the focus are shown in Fig. 3, highlighting how the algorithm predicts an error of up to 40 kPa near the focal region, which is very close to the Monte Carlo estimate.

A representative example of the empirical probability distribution for the pressure value at the focus is shown in Fig. 4. This shows a histogram of the pressure values at a single grid point for each of the 200 Monte Carlo simulations. The normal distribution predicted by the linear uncertainty propagation is shown on top. Despite the fact that the empirical distribution appears to be slightly skewed and, therefore, is deviating marginally from the supposed Gaussian distribution, the approximated distribution still satisfactorily explains the spread of possible pressure values for the given pixel.

Figure 5 shows the convergence of the estimated uncertainty using Monte Carlo for ten different random seeds. With small numbers of iterations (simulations), the Monte Carlo estimate has a large amount of variability. While the Monte Carlo solution will eventually converge to an unbiased value, this can take, on average, a very large number of simulations. In contrast, linear uncertainty propagation only requires two simulations, so it is computationally efficient. However, while for the simulation parameters used in this study it provides a good approximation of the true variance, for large values of *ε*, the biased nature of linear uncertainty propagation may start to become relevant.

## 5. Summary

We show that linear uncertainty propagation can be used to estimate uncertainty in simulated transcranial ultrasound fields when there is uncertainty in the medium property mapping between mass density and sound speed. Compared to Monte Carlo sampling, this provides a computationally efficient way of putting error bars on acoustic simulations used for treatment planning due to uncertainty in the material properties. In future, we will explore accounting for additional uncertainties (for example, in the skull attenuation) and the sensitivity of the predicted field to variations in different material properties and evaluate the regime of validity of linear uncertainty propagation for different model parameters (e.g., transducer shape, frequency, target region, etc.).

Another potentially interesting extension would be to directly couple the acoustic simulation with a thermal simulation to estimate the uncertainty on the temperature rise induced by the acoustic field (this is relevant for thermal ablation as well as safety for non-thermal therapies). As long as the thermal simulation is also written using a differentiable language, the algorithm outlined in this paper will still be applicable.

Last, note that the algorithm derived for estimating uncertainty is itself differentiable. This allows the predicted uncertainty map to be included in another optimization loop and the definition of a cost function that depends on uncertainty that is still optimizable using gradient descent, as done, for example, in Gerlach (2020) using a different method for estimating uncertainties.

## Acknowledgments

The authors would like to thank Taylor Webb and Kim Butts Pauly for provision of the raw data from Webb (2018) and Webb (2021) used to generate Fig. 1(b). This work was supported by the Engineering and Physical Sciences Research Council (EPSRC), UK, Grant Nos. EP/S026371/1 and EP/T022280/1.

## REFERENCES

*in situ*acoustic intensity estimates using MR acoustic radiation force imaging in combination with multifrequency MR elastography

*in vitro*results

*ε*-free inference of simulation models with Bayesian conditional density estimation