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.

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).

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.

Fig. 1.

(a) Conversion between Hounsfield units and mass density obtained using an image of a CIRS electron density phantom (model 062M). The error bars show the standard deviation in the image values averaged over each tissue equivalent electron density plug. (b) Conversion between density and sound speed using the data from Webb (2018) and Webb (2021). Each data point shows the density and sound speed measured for a specific skull sample (one outlier was removed from the data). The dashed line shows a linear fit to the data, and the solid lines show the 95% confidence interval for the linear model. The cross-hairs show the mean values.

Fig. 1.

(a) Conversion between Hounsfield units and mass density obtained using an image of a CIRS electron density phantom (model 062M). The error bars show the standard deviation in the image values averaged over each tissue equivalent electron density plug. (b) Conversion between density and sound speed using the data from Webb (2018) and Webb (2021). Each data point shows the density and sound speed measured for a specific skull sample (one outlier was removed from the data). The dashed line shows a linear fit to the data, and the solid lines show the 95% confidence interval for the linear model. The cross-hairs show the mean values.

Close modal
Using the data shown in Fig. 1(b), it is assumed the mapping between density ρ and sound speed c is given by the linear relationship
c ( α , β ) = c ̂ + α ( ρ ρ ̂ ) + β ,
(1)
where c ̂ = 2481 m·s–1 and ρ ̂ = 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 α 0 = 1.333 and β 0 = 0 m·s–1, with a standard deviation of σ α = 0.168 and σ β = 18.8 m·s–1.
Rather than performing a single acoustic simulation with a fixed mapping between density and sound speed (i.e., just using the values of α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 = ( α , β ),
p = f ( v ) .
(2)
Here, 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 μv with variance ∑v, given by
μ v = ( α 0 , β 0 ) , Σ v = ( σ α 2 0 0 σ β 2 ) ,
(3)
that is, v = μ v + ε, with ε N ( 0 , v ). If ε is small enough, Eq. (2) can be expanded as a Taylor series up to second order,
f ( v ) f ( μ v ) + J ( v μ v ) , with J = f | v = μ v = ( f / α f / β ) ,
(4)
where J is the Jacobian.
The truncated first-order Taylor series expansion of f ( v ) is a linear function of v; thus, inputs drawn from a normal distribution are mapped to outputs that are also normally distributed. The mean value of the distribution is f ( μ v ) (in other words, the simulated pressure field using α0 and β0), while the variance is given by Giordano (2016),
Var ( f ( v ) ) = J Σ v J T .
(5)
Since the variables are uncorrelated, the formula for the variance simplifies to
Var ( f ( v ) ) = ( f α | v = μ v ) 2 σ α 2 + ( f β | v = μ v ) 2 σ β 2 .
(6)
If the Jacobian is non-singular, this estimate can provide a good approximation for the variance for small variations of the parameter vector v. To estimate the normal distribution of the outputs, all that is needed is an efficient way to calculate f ( μ v ) and the Jacobian vector of f at v = μ v.

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 mm3, and truncated to a 120 × 70 × 70 mm3 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 f / α and f / β, 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.

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.

Fig. 2.

(a) The simulated pressure field, (b) the standard deviation of the simulation results using 200 Monte Carlo samples, (c) the standard deviation estimated by the linear uncertainty propagation method, and (d) the error between the estimated standard deviations. STD, standard deviation.

Fig. 2.

(a) The simulated pressure field, (b) the standard deviation of the simulation results using 200 Monte Carlo samples, (c) the standard deviation estimated by the linear uncertainty propagation method, and (d) the error between the estimated standard deviations. STD, standard deviation.

Close modal
Fig. 3.

The standard deviations estimated using Monte Carlo and linear uncertainty propagation, and the difference error between the two, for the (a) axial and (b) lateral lines passing through the focus. Subplots (c) and (d) show the axial and lateral pressure profiles through the focus, along with the 95% confidence interval (2-sigma) estimated using linear uncertainty propagation. STD, standard deviation; MC, Monte Carlo; Linear, linear uncertainty propagation.

Fig. 3.

The standard deviations estimated using Monte Carlo and linear uncertainty propagation, and the difference error between the two, for the (a) axial and (b) lateral lines passing through the focus. Subplots (c) and (d) show the axial and lateral pressure profiles through the focus, along with the 95% confidence interval (2-sigma) estimated using linear uncertainty propagation. STD, standard deviation; MC, Monte Carlo; Linear, linear uncertainty propagation.

Close modal

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.

Fig. 4.

The histogram of the empirical distribution of the pressure value at the focus, estimated using 200 Monte Carlo runs, and the corresponding normal distribution estimated by linear uncertainty propagation. LUP, linear uncertainty propagation.

Fig. 4.

The histogram of the empirical distribution of the pressure value at the focus, estimated using 200 Monte Carlo runs, and the corresponding normal distribution estimated by linear uncertainty propagation. LUP, linear uncertainty propagation.

Close modal

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.

Fig. 5.

Example of ten different Monte Carlo runs and their convergence in estimating the standard deviation at the focus, along with the uncertainty estimated using linear uncertainty propagation. MC, Monte Carlo; LUP, linear uncertainty propagation.

Fig. 5.

Example of ten different Monte Carlo runs and their convergence in estimating the standard deviation at the focus, along with the uncertainty estimated using linear uncertainty propagation. MC, Monte Carlo; LUP, linear uncertainty propagation.

Close modal

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.

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.

1.
Aubry
,
J.-F.
,
Bates
,
O.
,
Boehm
,
C.
,
Butts Pauly
,
K.
,
Christensen
,
D.
,
Cueto
,
C.
,
Gélat
,
P.
,
Guasch
,
L.
,
Jaros
,
J.
,
Jing
,
Y.
,
Jones
,
R.
,
Li
,
N.
,
Marty
,
P.
,
Montanaro
,
H.
,
Neufeld
,
E.
,
Pichardo
,
S.
,
Pinton
,
G.
,
Pulkkinen
,
A.
,
Stanziola
,
A.
,
Thielscher
,
A.
,
Treeby
,
B.
, and
van 't Wout
,
E.
(
2022
). “
Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models
,”
J. Acoust. Soc. Am.
152
(
2
),
1003
1019
.
2.
Caballero-Insaurriaga
,
J.
,
Rodríguez-Rojas
,
R.
,
Martínez-Fernández
,
R.
,
Del-Alamo
,
M.
,
Díaz-Jiménez
,
L.
,
Ávila
,
M.
,
Martínez-Rodrigo
,
M.
,
García-Polo
,
P.
, and
Pineda-Pardo
,
J. A.
(
2019
). “
Zero TE MRI applications to transcranial MR-guided focused ultrasound: Patient screening and treatment efficiency estimation
,”
J. Magn. Reson. Imaging
50
(
5
),
1583
1592
.
3.
Deffieux
,
T.
, and
Konofagou
,
E. E.
(
2010
). “
Numerical study of a simple transcranial focused ultrasound system applied to blood-brain barrier opening
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
57
(
12
),
2637
2653
.
4.
Gerlach
,
A. R.
,
Leonard
,
A.
,
Rogers
,
J.
, and
Rackauckas
,
C.
(
2020
). “
The Koopman expectation: An operator theoretic method for efficient analysis and optimization of uncertain hybrid dynamical systems
,” arXiv:2008.08737.
5.
Giordano
,
M.
(
2016
). “
Uncertainty propagation with functionally correlated quantities
,” arXiv:1610.08716.
6.
Hughes
,
T. W.
,
Williamson
,
I. A.
,
Minkov
,
M.
, and
Fan
,
S.
(
2019
). “
Forward-mode differentiation of Maxwell's equations
,”
ACS Photonics
6
(
11
),
3010
3016
.
7.
Lee
,
W.
,
Kim
,
H.-C.
,
Jung
,
Y.
,
Chung
,
Y. A.
,
Song
,
I.-U.
,
Lee
,
J.-H.
, and
Yoo
,
S.-S.
(
2016
). “
Transcranial focused ultrasound stimulation of human primary visual cortex
,”
Sci. Rep.
6
(
1
),
34026
.
8.
Li
,
N.
,
Gaur
,
P.
,
Quah
,
K.
, and
Butts Pauly
,
K.
(
2022
). “
Improving in situ acoustic intensity estimates using MR acoustic radiation force imaging in combination with multifrequency MR elastography
,”
Magn. Reson. Med.
88
(
4
),
1673
1689
.
9.
Marquet
,
F.
,
Pernot
,
M.
,
Aubry
,
J.-F.
,
Montaldo
,
G.
,
Marsac
,
L.
,
Tanter
,
M.
, and
Fink
,
M.
(
2009
). “
Non-invasive transcranial ultrasound therapy based on a 3D CT scan: Protocol validation and in vitro results
,”
Phys. Med. Biol.
54
(
9
),
2597
2613
.
10.
McDannold
,
N.
,
White
,
P. J.
, and
Cosgrove
,
R.
(
2019
). “
Elementwise approach for simulating transcranial MRI-guided focused ultrasound thermal ablation
,”
Phys. Rev. Res.
1
(
3
),
033205
.
11.
Miscouridou
,
M.
,
Pineda-Pardo
,
J. A.
,
Stagg
,
C. J.
,
Treeby
,
B. E.
, and
Stanziola
,
A.
(
2022
). “
Classical and learned MR to pseudo-CT mappings for accurate transcranial ultrasound simulation
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
69
(
10
),
2896
2905
.
12.
Montanaro
,
H.
,
Pasquinelli
,
C.
,
Lee
,
H. J.
,
Kim
,
H.
,
Siebner
,
H. R.
,
Kuster
,
N.
,
Thielscher
,
A.
, and
Neufeld
,
E.
(
2021
). “
The impact of CT image parameters and skull heterogeneity modeling on the accuracy of transcranial focused ultrasound simulations
,”
J. Neural Eng.
18
(
4
),
046041
.
13.
Papamakarios
,
G.
, and
Murray
,
I.
(
2016
). “
Fast ε-free inference of simulation models with Bayesian conditional density estimation
,” arXiv:1605.06376.
14.
Robertson
,
J.
,
Martin
,
E.
,
Cox
,
B.
, and
Treeby
,
B. E.
(
2017
). “
Sensitivity of simulated transcranial ultrasound fields to acoustic medium property maps
,”
Phys. Med. Biol.
62
(
7
),
2559
2580
.
15.
Stanziola
,
A.
,
Arridge
,
S. R.
,
Cox
,
B. T.
, and
Treeby
,
B. E.
(
2021
). “
A research framework for writing differentiable PDE discretizations in JAX
,” arXiv:2111.05218.
16.
Stanziola
,
A.
,
Arridge
,
S. R.
,
Cox
,
B. T.
, and
Treeby
,
B. E.
(
2022
). “
j-Wave: An open-source differentiable wave simulator
,” arXiv:2207.01499.
17.
Vaughan
,
T. E.
, and
Hynynen
,
K.
(
2002
). “
Effects of parameter errors in the simulation of transcranial focused ultrasound
,”
Phys. Med. Biol.
47
(
1
),
37
45
.
18.
Webb
,
T. D.
,
Leung
,
S. A.
,
Ghanouni
,
P.
,
Dahl
,
J. J.
,
Pelc
,
N. J.
, and
Pauly
,
K. B.
(
2021
). “
Acoustic attenuation: Multifrequency measurement and relationship to CT and MR imaging
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
68
(
5
),
1532
1545
.
19.
Webb
,
T. D.
,
Leung
,
S. A.
,
Rosenberg
,
J.
,
Ghanouni
,
P.
,
Dahl
,
J. J.
,
Pelc
,
N. J.
, and
Pauly
,
K. B.
(
2018
). “
Measurements of the relationship between CT Hounsfield units and acoustic velocity and how it changes with photon energy and reconstruction method
,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control
65
(
7
),
1111
1124
.
20.
Wise
,
E. S.
,
Cox
,
B.
,
Jaros
,
J.
, and
Treeby
,
B. E.
(
2019
). “
Representing arbitrary acoustic source and sensor distributions in fourier collocation methods
,”
J. Acoust. Soc. Am.
146
(
1
),
278
288
.