Elementwise approach for simulating transcranial MRI-guided focused ultrasound thermal ablation

This work explored an elementwise approach to model transcranial MRI-guided focused ultrasound (TcMRgFUS) thermal ablation, a noninvasive approach to neurosurgery. Each element of the phased array transducer was simulated individually and could be simultaneously loaded into computer memory, allowing for rapid (~2.5 s) calculation of the pressure field for different phase offsets used for beam steering and aberration correction. We simulated the pressure distribution for 431 sonications in 32 patients, applied the phase and magnitude values used during treatment, and estimated the resulting temperature rise. We systematically varied the relationship between CT (computerized tomography)-derived skull density and the acoustic attenuation and sound speed to obtain the best agreement between the predictions and MR temperature imaging (MRTI). The optimization was validated with simulations of 396 sonications from 40 additional treatments. After optimization, the predicted and measured heating agreed well (R2: 0.74 patients 1–32; 0.71 patients 33–72). The dimensions and obliquity of the heating in the simulated temperature maps were correlated with the MRTI (R2: 0.62, 0.74, respectively), but the measured heating was more spatially diffuse. The energy needed to achieve ablation varied by an order of magnitude (3.3–36.1 kJ). While this elementwise approach required more computation time up front (the combined simulation matrices were approximately 4.6 times higher than a single large simulation), it could be performed in parallel on a computing cluster. It allows for rapid calculation of the three-dimensional heating at the focus for different phase and magnitude values on the array. We also show how this approach can be used to optimize the relationship between CT-derived skull density and acoustic properties. While the relationships found here need further validation in a larger patient population, these results demonstrate the promise of this approach to model TcMRgFUS.


I. INTRODUCTION
Transcranial MRI-guided focused ultrasound has emerged as a noninvasive neurosurgical approach for thermal lesioning in the brain. This method uses a hemispherical phased array transducer to correct for aberrations of the acoustic field caused by the skull, allowing for accurate focusing to central regions in the brain. The method is clinically approved in several countries for thalamotomy for essential tremor [1] and tremor-dominant Parkinson's disease [2], and is under investigation for a number of other functional neurosurgery applications, including thalamotomy for neuropathic pain [3], pallidotomy for Parkinson's disease [4], and capsulotomy for obsessive-compulsive disorder [5]. Similar systems are being tested clinically to disrupt the blood-brain barrier as a treatment for Alzheimer's disease [6] and to enhance delivery of chemotherapy to patients with glioblastoma [7].
TcMRgFUS thermal ablation has several significant limitations. First, the treatment is currently restricted to central locations in the brain. When the transducer is focused at more peripheral targets, less energy can be transmitted through the bone due to large incidence angles between the transducer elements and the skull, leading to overheating. Furthermore, even at central locations, the shape of the focus can become oblique [8], which can put nearby structures at risk for unwanted thermal damage and side effects [9]. The bony properties of the skull vary substantially between patients, and the energy needed to achieve an ablative thermal dose at the focus can be excessive in some patients. As a result, some patients are not candidates for this treatment [10]. Finally, for reasons not fully understood, as sonications are delivered at escalating acoustic energies, it becomes more difficult to heat the tissue at the focus [11]. Being able to improve the focusing through the skull could ameliorate some of the challenges, expand the patient population who are candidates for TcMRgFUS and, perhaps, prevent side effects due to off-target heating.
To correct for aberrations induced by the irregular skull, the TcMRgFUS system estimates changes in of the ultrasound field induced by the bone for each element of the transducer. This correction is based on an acoustic model that uses the geometry of the skull and estimates of bone density obtained from CT (computerized tomography) scans [12,13]. The density is used to predict the acoustic sound speed and attenuation based on calibrations obtained with cadaver skulls [14][15][16]. Currently, a simplified proprietary acoustic model is used to generate the phase aberration corrections. Several studies have evaluated threedimensional acoustic models to predict the focal heating [8,[16][17][18]. Those studies have been able to reasonably predict the temperature and shape of the focus and the acoustic energy needed to reach an ablative thermal dose at the focus.
This work took a different approach to the simulation problem. Instead of simulating the entire acoustic field in one large model, we simulated each transducer element separately (Fig. 1). These simulations were rotated and interpolated into a global frame of reference and saved [19]. All the elementwise simulations could be loaded into computer memory simultaneously, enabling rapid (~2.5 s) iteration of phase and magnitude elements to simulate different correction schemes (or no correction) and perform beam steering to different targets. Ultimately, we aim to compile a library of these small simulations to produce a look-up table and to provide inputs for machine learning, which could enable us to rapidly simulate the heating during a treatment by relating each element to a previously simulated one. Rapid optimization of the phase and magnitude could maximize the peak intensity at the focus and better define the shape and size of the focal ablative region.

A. TcMRgFUS treatments
The treatments used the ExAblate Neuro transcranial MRI-guided focused ultrasound (TcMRgFUS) system (InSightec, Haifa, Israel), which operated at 660 kHz. This device has a 1024-element hemispherical phased array (993 active elements at our site) integrated with a 3T MRI (GE750, GE Healthcare). The patient was placed in a stereotactic frame which was attached to the MRI table. A flexible membrane was attached to the patient's head and the open face of the transducer. Acoustic coupling was achieved with degassed water that is chilled and circulated between sonications to minimize skull heating. The transducer was attached to a manually operated positioning system that was steered so that the focus was within ±1 mm of the initial brain target identified by the neurosurgeon. Additional steering from this location was achieved electronically using the phased array. The location of the transducer was found in the MRI space via MRI tracking coils, and imaging was performed using the body coil.
Before treatment, a CT scan of the head was obtained with a bone reconstruction kernel and a slice thickness of 1 mm or less. On the day of treatment, anatomic MRI scans were acquired in three orientations using a 3D FIESTA sequence [repetition time (TR)/echo time (TE): 5.1/2.4 ms; flip angle: 55°; receiver bandwidth: ±27.3 kHz; field of view: 22 cm; slice thickness 2 mm; matrix (typical) 224 × 288 × 28]. The CT and MRI were registered to each other using the software of the ExAblate system. Phase aberrations were calculated using a proprietary method by the manufacturer. The amplitudes of each phased array element were set by the device software to normalize the acoustic exposure over the head to avoid hotspots. The transformation between the CT and MRI space, the target locations in MRI space, and the magnitude and phase of each transducer element were saved for each treatment. This information was used to simulate the acoustic field as described below.
During each sonication, magnetic resonance temperature imaging (MRTI) was obtained in a time series in a single, operator-defined plane using a fast spoiled gradient echo sequence. The orientation of the imaging plane and frequency/phase encoding directions varied between each sonication. Two acquisitions were obtained before the start of each sonication; the first was ignored due to artifacts. Acquisition continued over the course of the sonication; five additional acquisitions were obtained to map cooling. In the first 17 patients, MRTI parameters were TR/TE: 27.8/12.9 ms; flip angle: 30°; receiver bandwidth ±5.7 kHz; field of view: 28 cm; slice thickness: 3 mm; matrix: 256 × 128. In the others, a multiecho readout was performed (TR: 28.0 ms; TE: 3.1/7.8/12.5/17.2/21.9 ms; flip angle: 30°; receiver bandwidth: ±35.71 kHz; field of view: 28 cm; slice thickness: 3 mm; matrix: 256 × 128). The higher readout bandwidth of the multiecho sequence reduced spatial distortions in the frequency-encoding direction with minimal loss of SNR [20].
The MRI reconstructed magnitude, real and imaginary images that were converted to phasedifference images. Temperature changes were estimated from the phase-differences [21] using a temperature dependence of −0.00909 ppm/°C [22]. This value was selected because it is what is used in the device software. The raw images were stored in DICOM format and independently analyzed offline using MATLAB software developed in-house. Artifacts resulting from changes in magnetic field (presumably due to patient motion outside the brain) were removed using a procedure outlined elsewhere where nonheated regions in each phase-difference images were fit to a smooth 2D surface [23]. The surface was extrapolated into the heated region and subtracted off. The phase difference maps obtained with the multiecho sequence were combined via a weighted average based on the expected signal-tonoise ratio of the MRTI [24], with the weight for the nth echo given by W n = TE n e − TE n T 2 * 2 . (1) We used 30 ms for the T 2 * relaxation time for brain.

B. Simulations
Numerical modelling of the acoustic field was performed using k-Wave, an open-source MATLAB toolbox [25]. For every patient, the pressure amplitude was simulated individually for each element of the array. The element was assumed to be a 1 cm diameter flat circular piston centered on the location provided by the manufacturer. When all simulations were loaded into memory, applying magnitude and phase and calculating the total pressure amplitude took approximately 2.5 s. Phase corrections for "ideal" focusing were found by simply subtracting the phase at the target location from each element's simulated pressure.
Next, we used the bioheat equation to estimate the heating [28]. We averaged the heating to match the voxel dimensions and temporal resolution of the MRTI with the imaging orientation used during treatment.
The parameters used in the numerical modelling of the acoustic pressure and the subsequent temperature rise are listed in Table I. To estimate the sound speed and attenuation of the skull, we first used the density estimated from the CT scans and the empirical relationships presented by Pichardo et al. [15] interpolated to 660 kHz. We assumed a linear relationship between Hounsfield units and density, and −1000 and 57 Hounsfield units for air and soft tissue, respectively.
Overall, we considered 72 patient treatments. The first 32 patients were examined in detail and were used to optimize the relationships between skull density and the acoustic attenuation and sound speed. Thirty of these patients were treated in the ventral intermediate (VIM) nucleus of the thalamus for essential tremor; two were treated in the globus pallidus for Parkinson's disease. The number of sonications per treatment and the acoustic parameters used (power, duration) varied between patients. Details for each patient, their treatments, and the CT scans are listed in Table II. A total of 447 sonications were performed in the treatments in patients 1-32; 16 sonications were excluded due to artifacts in the focal region in MRTI. After completing the analyses of the first 32 patients, we validated our results by simulating 40 additional essential tremor patients (Table S1).
The simulations were performed in MATLAB using the O2 High Performance Compute Cluster, supported by the Research Computing Group at Harvard Medical School. We stored the complex steady-state pressure amplitude after each simulation. The simulations in the first 32 patients were run five times: using the density/attenuation relationship described in Pichardo et al. [15], without attenuation, using the optimized density/attenuation and finally with the two optimized density/sound speed relationships. The optimization procedure is described below. We simulated the additional 40 patients with the two optimized density/ sound speed relationships and without attenuation. The elemental simulations typically required approximately 20 min each to run.

C. Optimizing skull attenuation
We investigated the feasibility of systematically iterating the elemental simulations to find a relationship between attenuation and skull density that resulted in heating that better matched the MRTI for the first 32 patients. It was impractical to repeatedly run the full numerical model, so we used a simplified attenuation model that we applied to simulations performed without attenuation. This simplified model assumed that the primary attenuation effects occur along the z direction: where P k0 (x i , y i , z i ) is the pressure distribution simulated without attenuation for element k, α(x i , y i , z i ) is the attenuation, and Δz is the grid size.
We assumed the relationship between attenuation and skull density p sk (x i , y i , z i ) can be approximated by a series of polynomials: resulting in where N b (x i , y i , z i ) is the number of soft tissue voxels between the transducer and z i ,. After exchanging the order of summation this becomes The values for N b (x i , y i , z i ) and the cumulative summations of the skull density (∑ j = 0 i ρ sk x i , y i , z j m ) for orders m = 0 -4 were calculated for each transducer element, interpolated to the array space, and saved. These data, along with the individual pressure estimates for each element could then be loaded into memory. The total three-dimensional pressure amplitude for all elements and the resulting temperature rise could then be estimated for different values of A m . We only included the 23 × 23 × 31 spatial simulation points centered on the focus location of the simulation to save time.
We attempted to find the A m coefficients that minimized the difference between the simulated and measured peak temperature rise for all patients. It was not practical to load the simulated pressure fields and skull density summations for the 32 patients into computer memory simultaneously. We thus systematically explored different coefficients for each patient separately and in parallel. We defined a grid of nodes at 12 values of density between 1200 and 3500 kg/m 3 and 10 values of attenuation between 80 and 500 dB/cm (Fig. 2). We then selected a set of three nodes, which were each then fit to a fourth-order polynomial (order chosen ad hoc) to determine the A m coefficients. Two of these nodes had densities of 1200 and 3500 kg/m 3 (100 possible sets), and the third was one of the 100 nodes between these two. Thus, for each patient we examined 10 000 sets of A m coefficients.
For each set of coefficients, we used Eq. (5) to estimate the attenuated pressure field for the individual elements. We then set the phase and magnitude for each sonication and estimated the focal temperature rise. After compiling these temperature estimates, we found the coefficients that minimized the mean squared error between the simulated and measured temperature rise. Finally, we repeated the full simulations with the optimized attenuation/ density relationship.
As shown in the results below, in many patients, plots of the difference between measured and simulated heating as a function of previously applied energy was relatively constant for the initial sonications, after which it began to deviate. We interpreted this deviation as an irreversible change in skull properties. Thus, when we found the optimized A m coefficients, we only included sonications that occurred before this deviation.

D. Optimizing skull sound speed
The deviations between simulations and MRTI in plots of heating as a function of acoustic energy suggest that the skull acoustic properties might be changed during treatment.
Assuming ad hoc that this change was primarily in sound speed, we investigated whether we could find optimized relationships between sound speed and skull density that could be used before and after this change occurred. As described above for attenuation, we created a simplified model that allowed us to rapidly explore different density/sound speed relationships using previously obtained simulations. We assumed that the relationship between the inverse of the sound speed and the skull density can be approximated by a series of polynomials: 1 c sk x i , y i , z i = ∑ m B m ρ sk x i , y i , z i m . (6) To estimate the effect of changing the skull sound speed on the transmitted pressure magnitude, we assumed the biggest loss occurred due to the initial reflections at the interfaces between soft tissue, the inner and outer tables, and the diploe (i.e., we ignored multiple reflections and other reflections within the skull). The transmission coefficient was calculated using the following relationship [12]: where the acoustic impedance of the outer and inner tables, the diploe, and soft tissue, respectively, are given by We used brain values for the soft tissue density (ρ T ) and sound speed (c T ), and density values measured at the outer and inner tables of the skull and the diploe for ρ OT (x i , y i ), ρ IT (x i , y i ), and ρ DP (x i , y i ), respectively, as described below. The incident and transmitted angles with respect to the normal vectors at the outer and inner tables and the diploe (θ i , θ t , θ i ′, θ t ′, respectively) were found via Snell's law: where γ O , γ dp , and γ I were the incidence angles of the outer skull surface, the center of the diploe, and the inner skull surface measured with respect to the z axis. A diagram showing how the incidence angles were defined for this procedure is shown in Fig. S1. We measured the density and angles of incidence for the inner and outer tables with respect to the z axis, for each (x i , y i ) coordinate as described below.
To estimate phase shifts resulting from changing the skull sound speed, we multiplied the delays expected the along the z direction by 2π times the FUS frequency (f): After combining with Eq. (6) and changing the order of summation, this becomes The cumulative summations of the skull density ∑ j = 0 i ρ skull x i , y i , z j m saved for the attenuation estimation were used.
We applied the magnitude and phase changes to simulations performed with the optimized attenuation. We used Eqs. (7) and (11) to remove the effects of the density/sound speed relationship used in those simulations before testing new ones.
As was done with the attenuation optimization, we estimated the temperature rise over a predefined set of density/sound speed curves for every sonication. We defined a grid of nodes covering densities between 1200 and 3500 kg/m 3 and sound speeds between 1500 and 4000 m/s. We then selected a set of three nodes, and we fit the inverse of the sound speed to a fourth-order polynomial to determine the B m coefficients. Two of these nodes had densities of 1200 and 3500 kg/m 3 (100 possible sets), and the third was one of the 100 nodes between these two. Thus, for each patient, we examined 10 000 sets of B m coefficients. The relationship that minimized the error between the simulations and measurements was found for the sonications that deviated from the model after attenuation optimization. We then repeated the full simulations with the optimized attenuation and sound speed relationships.

E. Data analysis
All analysis was performed in MATLAB. We compared the peak temperature rise at the focus measured with MRTI to that estimated with the simulations using linear regression. To compare "treatment efficiency", we estimated the energy needed in each patient to achieve a focal temperature of 55 °C. We noted that plots of measured focal heating as a function of acoustic energy generally followed a logarithmic curve. Thu, we used nonlinear least squares regression (using the function "nlinfit" in MATLAB) to estimate this energy. To further test our ability to predict the "treatment efficacy", we compared the measured and predicted slopes of plots of the logarithm of the acoustic energy and the focal temperature rise.
We also compared the shape and orientation of the heating. The simulated and measured temperature maps at the time of peak temperature rise were fit using the Levenberg-Marquardt method via the MATLAB function "fit.m" to a two-dimensional Gaussian distribution: ΔT (x, y) = A + B exp − x − x 0 cosθ + y − y 0 sinθ σ x 2 × exp − x − x 0 sinθ + y − y 0 cosθ σ y 2 . (12) We used linear regression to compare the resulting estimates of the dimensions (σ x , σ y ) and obliquity (θ). For the obliquity, not all the temperature maps had a well-defined angle. Thus values for θ that had 95% confidence intervals greater than ±30° were excluded. In each patient, we selected representative sonications that were largely free of artifacts in the MRTI.
We selected examples from each imaging orientation and frequency encoding direction for each patient. Not all orientations and encoding directions were available in every patient. Overall, we examined 110 sonications when comparing the focal heating shape.
We examined whether different factors gleaned from the skull CT could predict the acoustic energy needed to achieve an effective thermal exposure at the focus (i.e., energy required to heat to 55 °C). These factors were determined with the individual volumes used in the elementwise simulations (Fig. 3). Working along the direction of propagation, we identified the first and last skull voxel, the voxels at the center of the inner and outer tables, and the voxel that had the minimum density between the outer and inner table for each of the 44 × 44 coordinates in the elementwise simulation. The coordinates of the first and last skull voxel and the diploe were fit to smooth surfaces. Normal vectors to this surface were calculated to obtain the angles of incidence of the skull and the direction of ultrasound propagation. We used these data to calculate different skull-derived factors (see Ref. [29] for more details). We then examined the relationship between each factor and the energy needed to reach a peak focal temperature of 55 °C. The ability of different metrics to predict this energy were evaluated using linear regression with the "fitlm.m" function in MATLAB.
Since there were obvious outliers, we used the robust fitting option which uses iteratively reweighted least squares with a bisquare weighting function.

III. RESULTS
A. Comparing shape, size, and obliquity of simulated and measured temperature maps This general agreement was evident overall. Figure 5(a) shows example simulated temperature maps and corresponding MRTI for patients 1-32 in different imaging orientations. The relative size and obliquity of the focal heating were similar in most patients. However, the heating was more diffuse in the MRTI than was predicted by the simulations; Fig. 5(a) shows 50% isotherms in the MRTI and 25% isotherms in the simulations.
To characterize the ability of the simulations to predict the shape of the heating measured with MRTI, we estimated the dimensions and obliquity by fitting the focal heating in the MRTI and the simulated temperature maps to a two-dimensional Gaussian distribution [Figs.
. A correlation was observed between both the dimensions and tilt angles measured in the MRTI correlated and those predicted by the simulations [correlation coefficient (R 2 ): 0.62 and 0.74, respectively; probability value (P < 0.001)]. However, the dimensions of the simulations were smaller than the measurements in most cases. The simulations predicted a tilt in the same direction as the MRTI, but the absolute angle was less in the simulations, particularly in the left/right direction in coronal MRTI. The ratios of the dimensions (σ x /σ y ) of these fits for the simulations and measurements were correlated (R 2 : 0.76).
The simulated heating shown in Fig. 5 used the optimized relationships between density and acoustic properties described below. Similar results were found using the relationships described by Pichardo et al. [15]. There was no significant difference (P > 0.05) between the optimized and Pichardo relationships in the resulting dimensions or angularity of the simulated focal region.
Left/right obliquity and elongation in the superior/inferior direction poses a risk during thalamotomy of damaging the internal capsule, which is located laterally and inferior to the location of the thalamic target. An extreme example of such a case is shown in Figs. 6(a)-6(h). In this patient, there was obliquity in both the left/right and superior/inferior directions, leading to heating in the internal capsule that could not have been detected in the singleplane MRTI in any orientation [Figs. 6(c) and 6(f)]. This double obliquity was predicted by the simulations [Figs. 6(g)-6(h)]. Another example of this risk is shown in Figs. 6(i) and 6(j). The device manufacturer added an option to change the magnitude distribution of the phased array to reduce the left/right obliquity. Figures 6(i) and 6(j) show MRTI acquired with and without this magnitude mask. While a corresponding reduction in obliquity is also evident in the simulations, the shape of the predicted heating is less diffuse and does not include the lateral heat spread into in the internal capsule [arrow in Fig. 6(i)].

B. Predicting focal heating
The examples presented above demonstrate that the simulations generally predicted the relative size and shape of the focal regions that were observed in MRTI. We also compared measurements and simulations of the peak temperature at the focus for the 431 sonications delivered without MRTI artifacts in these 32 patients.
Our first simulations used the relationship between density and attenuation found by Pichardo et al. [15]. After comparing simulations and measurements (Fig. S2), we made two general observations. First, while in some cases good agreement was observed, we found that the simulated temperatures were generally less than the measurements, particularly when the exposure level was relatively low. Second, we noted that plots of acoustic energy versus heating followed parallel trajectories during the initial sonications performed during the treatments. However, in many patients, these trajectories deviated as the acoustic energy increased, and the measured heating increased less than the simulations predicted as the energy increased.
One interpretation of these two observations is that the relationship between skull attenuation and density was incorrect, and that in some patients the acoustic properties of the skull changed after a certain level of exposure, leading to defocusing. Thus we evaluated the feasibility of optimizing the density/attenuation relationship to better match the MRTI for those sonications delivered before this change. Figures 7(a) and 7(d) compare focal heating in MRTI to simulations using the attenuation/density relationship found by Pichardo et al. We next tested whether we could find a relationship between sound speed and density that predicted focal heating that better matched the measurements for the sonications both before and after the presumed change to the skull. The resulting estimated peak temperatures for the deviant sonications are shown with the red crosses in Figs. 7(c) and 7(f). The agreement with the measurements was improved on average, but the heating in some patients was still not well predicted. Figure 8(a) shows plots of the measured and simulated temperature rise as a function of the acoustic energy for patients 1-32. With the optimized attenuation/density relationship, excluding the deviants, the trend in these plots was similar for measurements and simulations, and the simulations captured effects such as beam steering, changing sonication duration, and varying the plane used for the MRTI. These effects are evident in Fig. 8(a) where the simulated temperature/energy relationships are not smooth. The simulations grossly underpredicted the focal heating in a few patients (such as patients 2, 3, and 16). Furthermore, in most patients the measurements followed a logarithmic curve [appears linear in Fig. 8(a), which has a logarithmic scale on the x axis], while the simulations predicted a more linear relationship [appears curved in Fig. 8(a)]. Patient 6 is an example of this behavior. Figure 9 compares the slopes of the measured and predicted heating versus acoustic energy curves in Fig. 8(a). A good correlation between measurements and simulations was observed (R 2 : 0.57), but the simulated efficiency was slightly higher in most patients.
Using the optimized density/sound speed relationship improved the prediction of the focal heating for the deviant sonications. The predicted focal heating with this optimized relationship is shown in Fig. 8(a) by the green "×" symbols. There were some patients (patients 29 and 30, for example) that were still not well-predicted by the simulations after this optimization.
As we gained experience with these treatments, we used larger steps in acoustic energy between sonications. Thus, in two cases (patients 27 and 30), all sonications except the first few low-energy ones had heating less than the simulations predicted. In these two patients, there was a large difference between the simulation and the measurements even though the prior energy delivered was low.

C. Focal heating versus time
We compared plots of the temperature rise as a function of time for the measured and simulated MRTI [ Fig. 8(b)]. For each patient we selected two sonications, one where the heating was approximately 10 °C, and another where the temperature rise was approximately 18 °C (corresponding to an absolute temperature of 55 °C). The heating rates of the simulations and measurements were similar. However, the decay of the heating after the sonication was faster in the simulations than the measurements. To characterize this, we fit the temperature decay to an exponential function and estimated the time required to reduce the heating by one half. The mean half-life was 14.0 ± 3.7 s in the measurements, significantly higher (P < 0.001) than the 8.6 ± 1.1 s estimated in the simulations.

D. Optimizing relationships between skull density, attenuation, and sound speed
In optimizing the attenuation/density relationship, we used a simplified attenuation model that allowed us to rapidly estimate the focal heating for different scenarios. We estimated the temperature rise for every patient and sonication for 10 000 different attenuation versus density curves. The best result we found [excluding the deviants shown in Fig. 7(a) and Fig.  7(d)] is shown in Fig. 10(a). The relationships that were within 10% of the best curve followed a similar trajectory except for higher densities that were rarely observed in these 32 patients. Using the best relationship, the full simulations were repeated. Comparison of the focal heating between the simplified and full simulations are shown in Fig. 10(b). Overall the two estimates were correlated (R 2 : 0.94), with heating estimated by the simplified model slightly higher in most patients than that predicted by the full simulations.
After the attenuation optimization, we performed a similar procedure with sound speed. Figure 10(c) shows the density/sound speed relationships that best predicted the measured peak focal heating for the deviant and nondeviant sonications. Figure 10(d) compares the heating predicted using the simplified sound speed model to that predicted by the full simulation. The two were correlated (R 2 : 0.84), but there was considerable mismatch in some patients, and on average the heating was slightly higher in the simplified model.
The coefficients characterizing the optimized density/attenuation and sound density/speed found in Eqs. (3) and (6) that resulted in the best agreement between the models and the measurement are listed in Table III. The difference in heating between measurement and the simulation varied substantially for the different relationships tested (Fig. S3).

E. Validating the optimization
We simulated an additional 40 ET patients (396 sonications) using the optimized relationships we found between skull density and the attenuation and sound speed. Overall the agreement between simulation and measurement in the peak temperature rise was similar for these patients for the nondeviant sonications (Fig. 11). The optimization of the density/ sound speed relationship for the deviant sonications improved the prediction in many cases, but there were several cases where the simulations continued to overpredict the focal heating. Additional information about these patients is shown in Table S1 and Fig. S4. When the 72 patients were considered together (825 sonications), regression of the simulated versus measured temperature rise at the focus yielded a slope of 1.05 ± 0.02, and a y intercept of −1.31 ± 0.39 (R 2 : 0.71). The root mean square difference between simulation and measurement was 3.4 ± 3.1 °C.
We also repeated the optimization of the density/attenuation and density/sound speed relationships in these 40 patients and for all patients (Fig. 12). The best relationships found were similar to those from patients 1-32.

F. Effects of spatial averaging and aberration correction
We investigated the effects of spatial averaging over the imaging planes used for MRTI [ Fig.  13(a)] and aberration correction [ Fig. 13(b)] on the simulated temperature maps. Overall, the predicted temperature rise was 1.5 times higher before spatial averaging when the phase corrections used during treatment were used. Given that the measured heating was more diffuse than the simulations predicted, this may be overstating the effects of averaging. The simulations predicted that using the ideal phase corrections would result in improved heating 1.2 times higher than the corrections used during the treatments. The simulations also predicted that without the phase aberration corrections used during the treatments that the temperature rise would have been reduced by half on average.

G. Predicting ablative energy
There was substantial variability among the different patients in the ultrasound exposure levels needed to achieve a temperature rise sufficient to produce a thermal lesion, presumably due to differences in scattering or other differences in internal skull properties.
To characterize this variability, we estimated the acoustic energy required in each patient to achieve an ablative thermal exposure of 55 °C in MRTI [dotted vertical lines in Fig. 8(a)]. This energy range spanned an order of magnitude; it ranged 3.7-34.6 kJ in patients 1-32 and 3.3-36.1 kJ in patients 33-72. We compared the ability of the simulations to predict this exposure level along with multiple other metrics gleaned from the CT scans in patients 1-32 (Fig. S5). Correlation between the measured and simulated energy needed to reach 55 °C was significant (P < 0.001; R 2 : 0.78). Several other skull-derived metrics also showed a significant correlation. The three cases that required the highest energy were outliers in these regressions, and the high acoustic energy required to reach an ablative thermal exposure at the focus (>20 kJ) was only predicted by the simulations. Figure 14 shows color-coded maps of the skulls for patients 1-32 showing the spatial distribution of three metrics that had highly significant correlations with the energy needed to reach 55 °C: loss in pressure due to skull attenuation, skull thickness, and cortical/ trabecular density ratio. The patients are ordered by the acoustic energy needed to reach 55 °C. While trends are evident, there was substantial variability. Note, for example, patient 19, who's skull was relatively thick and highly attenuating but required very little energy, and patient 30, who's skull was thinner and less attenuating than others that required less energy.
We also compared these skull-derived metrics for patients that did and did not have sonications with heating/energy trajectories that deviated from the predictions. The skulls of the patients that had sonications with deviations were thinner (6.8 ± 1.0 versus 8.0 ± 1.1 mm; P < 0.001) and had a smaller volume (278 ± 47 versus 330 ± 62 cm 3 ; P < 0.001) than those where the trajectories did not deviate. The corresponding loss due to attenuation was also less (0.69 ± 0.06 versus 0.79 ± 0.74 ± 0.05; P < 0.001) in the patients with deviants. The standard deviation of the attenuation along the direction of ultrasound propagation was also significantly higher (67.5 ± 4.7 versus 60.7 ± 4.3 Np/m; P < 0.001).

H. Comparing "ideal" and treatment phase aberration corrections
The phases used for aberration correction determined by the device software and used in the treatments were grossly similar to the "ideal" phases predicted by the simulations in every case in patients 1-32 (R 2 : 0.84, P < 0.001; Fig. S6). We examined whether the different metrics gleaned from the CT scans were correlated with differences in the predicted ideal phase corrections and those used during treatment. Factors that had a significant correlation with the difference between the ideal and treatment phases were the number of elements with incidence angles greater than 25° (P < 0.001), the standard deviations of the density, sound speed, and impedance, the density and sound speed of the outer table (P < 0.05), and the difference between the internal and external incidence angles (P = 0.001). The proprietary model used by the manufacturer considers shear mode transmission for incidence angles above a certain threshold, so it is not surprising that differences that were correlated with factors related to incidence angles.

IV. DISCUSSION
This work describes our experience testing a simulation framework that enables rapid (~2.5 s) calculation of the pressure field for different transducer magnitude/phase distributions and that could be iterated using simplified models to identify relationships between the CTderived density and the skull acoustic properties that better matched the MRTI. This approach required more upfront work to generate the individual elemental simulations, but it could be performed in parallel using a computing cluster and potentially allows for more flexibility in exploring how changing different properties affects focal heating.
Ultimately, we aim to assemble a look-up table of elemental skull segments and pressure distributions to enable us to accurately predict the three-dimensional pressure distribution in real time. This table could be based on similarity comparisons between skull segments or based on derived factors such as thickness, angles, etc. Machine learning methods may be useful in developing these predictive models to relate the skull segments to the resulting pressure distributions. If one had the ability to predict the three-dimensional pressure distribution in real time, one could rapidly iterate the phase and magnitude of the transducer elements (or the position and angulation of the transducer itself) to optimize not only the peak temperature rise but also the shape and size of the focal heating. With an ability to predict the pressure amplitude in real time, one could also apply holographic methods [30] to shape the focal region to match the desired anatomy of interest. Such shaping may be more useful in TcMRgFUS applications other than thermal ablation, such as blood-brain barrier disruption [6,7], neuromodulation [31], and nonthermal ablation [32] that do not require high acoustic exposure levels, as modifying the spatial pressure field would likely reduce the peak pressure amplitude at the focus.
To achieve this goal, we need to better predict the shape of the focal heating. The heating measured with MRTI was more diffuse than the simulations predicted. This result could suggest that acoustic parameters used in the model were not correct, or it could reflect factors that were not considered here in the simulations. Transmission after shear mode conversion in elements with oblique incidence angles, for example, was not accounted for in these simulations. While the high shear attenuation coefficient will limit this transmission, studies in cadaver skulls suggest that transmission of up to 17%-23% of the power can be transmitted via a shear wave for incidence angles greater than 30° [33]. Complex reflections, such as those that occur from the face of the transducer and those that occurred outside of the elementwise simulation volumes may have also contributed, and probably to a small effect [34], nonlinear acoustic propagation. These components of the transmitted acoustic wave were not taken into account in the aberration correction, perhaps leading to diffuse heating at and around the focal region in the measurements that was not present in the simulations. This diffuse heating may explain the slower temperature decay after the sonications [ Fig. 8(b)]. A faster decay in the simulation might also suggest that the thermal properties and the perfusion coefficients used for brain were incorrect.
Other than the diffuse heating, the simulations did a good job overall in predicting the relative peak temperature rise and orientation/obliquity of the focal region. However, it is possible that the parameterization of the CT scans and the acoustic and thermal properties need refinement. For example, based on earlier works [14,15], we used a simple linear extrapolation to estimate skull density from Hounsfield units, and the Hounsfield units were estimated directly using information from the DICOM headers. Others have suggested potentially more accurate approaches to estimate skull density from CT scans [16,35], and have shown that scans from CT scanners from different vendors and with different reconstruction kernels can yield different results [36]. Here, the CT scans were from a variety of vendors, but we did not see any obvious effects on the simulations.
We began this project using the experimentally derived skull density/sound speed and density/attenuation relationships found by Pichardo et al. [15]. While we did observe good agreement between the measurements and the simulations in some patients, in most cases the peak temperature was underpredicted by the numerical model for sonications at relatively low exposure levels. Furthermore, in many patients, the relationship between acoustic energy and temperature rise followed a similar trajectory in the initial sonications but appeared to deviate after a certain amount of acoustic energy was applied.
Based on these results, as well as clinical experience showing that increasing acoustic energy often produces diminishing increases in heating as the treatment progresses [11] and that has observed skull damage in some patients after treatment at high energies [37], we hypothesized that (1) the relationship between attenuation and skull density found by Pichardo et al. was not optimal, and (2) the deviation could represent an irreversible change in skull acoustic properties. We attempted to find a density/attenuation relationship that better predicted the MRTI measurements for sonications delivered before the skull changes. This process improved the predictive ability of the simulations substantially, but there were still patients where the agreement was poor. It is possible that there is no single relationship between density and attenuation, as we do not capture the microstructure of the bone below the resolution of the CT scan, which likely affects the scattering of the acoustic field [38,39]. More work is needed to understand how well we can use clinical imaging to predict the attenuation for an individual patient.
We also explored whether we could find a skull density/sound speed relationship that could be used after a (presumed) change in skull properties. While we were able to identify a relationship that improved the predictive ability of the model after the simulated heating deviated from the measurements, there was still substantial variability. Validation of these relationships in 40 additional patients revealed similar agreement between predicted and measured heating overall. However, the density/sound speed relationship found in patients 1-32 after the presumed skull changes had even more variability in patients 33-72. Finding such variability is probably not surprising, since heat-induced changes to the skull will not be spatially uniform or simply binary. It may not even be correct to assume that the sound speed is the dominant factor that changes during the treatment.
Furthermore, our interpretation of the results that the skull acoustic properties changed may not be correct for every patient. Note, for example, patients 27 and 30 in Fig. 8(a), where the predicted heating with the optimized density/attenuation relationship was substantially higher than the measurements even though the previously delivered energy was low. Similar findings were observed in several treatments of the 40 additional patients (Fig. S4). These results might suggest that the deviations in temperature/energy were not due changes in skull properties from previously applied sonications, but perhaps instead a dynamic change in acoustic properties with heating. It could also be that the sonications in these patients should not have been included as "deviants" and that our optimized density/attenuation relationship needs further refinement. Since we do not have heating measurements at intermediate energies in these patients, we cannot answer this question. Future studies with additional patients are needed to understand cases like these.
Other factors may contribute as well. For example, the performance of the transducer elements might change if they heat at high power levels, or strong reflections from the skull could change their impedance.
We have demonstrated how one could use elementwise simulations to improve the relationships between acoustic attenuation and sound speed. Here, we used simplified models that could be rapidly applied to previously obtained simulations to explore the effects of thousands of different relationships on resulting focal heating. While these simplified models did a good job in predicting the temperature at the focus for the full simulation, a better approach would be to explore these relationships to match not only the peak temperature but also the shape and size of the focal heating. Perhaps it would also be better to optimize over multiple acoustic and thermal parameters simultaneously. Individual patients or groups of patients may have skulls with different relationships between density and acoustic properties. Repeating this analysis in more patients may reveal such trends. Finally, we need to evaluate whether the relationships between density and acoustic properties found here can be used prospectively improve the focusing in future TcMRgFUS treatments.

A. Predictors of treatment efficacy
We observed an order of magnitude range in energy required to achieve an effective level of heating. It would be useful to predict this energy before treatment to screen potential TcMRgFUS patients. Currently, the device manufacturer uses the "skull density ratio", which is the ratio in density between the diploic trabecular and dense cortical bone [10,40] to predict which patients will require high acoustic energies to achieve focal temperatures sufficient to induce a thermal lesion. We examined different skull-derived parameters to see how they compare in this prediction. The manufacturer-derived "skull density ratio" as well as our own calculation of this ratio were both predictive of the energy needed to reach 55 °C (R 2 : 0.35, 0.40, respectively; P < 0.01), but a number of other parameters, such as the skull thickness and the expected loss due to attenuation were also similarly predictive. However, there was significant variability and there were outliers that were not predicted by simple skull measurements. The simulations did a better job in predicting the acoustic energy needed to reach 55°C (R 2 : 0.78; P < 0.001) than any single factor.

B. Conclusions
This study demonstrates the feasibility of using an elementwise approach to simulate TcMRgFUS thermal ablation and predict the shape and magnitude of the focal heating in multiple orientations. While there was significant variability, the numeric model predicted the relative shape and focal temperature rise on average. Deviations between the measured and simulated heating were observed in many patients over the course of the treatments, perhaps reflecting a change in skull properties. We also demonstrated how this approach could be used to optimize the relationship between CT-derived density and the acoustic attenuation and sound speed. We found initial estimates of these relationships based on 32 patient treatments, including an estimated density/sound speed relationship that could be used after the model and the measurements deviated. The optimization was validated in 40 additional patients. Future work will expand on this study to refine these optimized relationships. Importantly, the optimized relationships need to be validated prospectively during clinical TcMRgFUS treatments.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. TcMRgFUS phased array is superimposed. A small volume is selected that includes one array element. Right: The pressure field (logarithmic scale shown) is simulated in this volume using the CT scan to estimate the skull acoustic properties. The simulation is then rotated and interpolated to a global volume that includes the focal region and saved. After repeating for every element, the combined field can be rapidly obtained with different phase/ magnitude values for the elements.

FIG. 2.
Optimizing the relationships between skull density and attenuation. We examined 10 000 relationships per sonication (black curves). To generate the curves, we used two nodes with densities of 1200 and 3500 kg/m 3 (red circles) and a third indicated by one of the yellow crosses. These nodes were fit to a polynomial. We generated a set of 100 curves for each yellow cross; an example set is shown in magenta.

FIG. 3.
Feature extraction from the CT scans for each phased array element. (a) For every x and y coordinates of the elemental simulation, threshold analysis of the density along the z direction was used to find the first and last voxel of the inner and outer skull surface. Working from the outer surface, the first peak encountered was used to identify the outer table. The same procedure was used working from the inner surface to find the inner table.
The minimum density between the inner and outer table was used to identify the diploe. The coordinates of the inner and outer surface and at the diploe [(b) and (c)] were fit to smooth surfaces (d). The normal vectors for these surfaces were found for each x and y coordinates (e); the individual normal vectors (black) and the mean vector (red lines) are shown. Note that for clarity, (b), (c), and (e) only plot every fourth x and y coordinates.

FIG. 4.
Example simulated pressure distributions (a) and measured and predicted temperature maps (b) from a pallidotomy treatment (patient 11). In (a), three orientations with respect to the transducer are shown for one sonication with no aberration correction, the correction used during treatment (Tx), and ideal correction; the simulated field without a skull (but with magnitude shading and beam steering) is also shown. The location of the geometric focus is indicated. For this sonication, the focus was electronically steered away from the geometric focus to x = 3 mm and y = 1 mm. In (b), MRTI for three sonications are shown; the insets are the simulated heating with the aberration correction used during treatment for these examples. Scale bars: 1 cm.

FIG. 5.
Comparing predicted and measured focal shape, size, and orientation. (a) Simulated temperature maps (Sim) and corresponding MRTI acquired during TcMRgFUS in 32 patients. One example from each imaging orientation was selected per patient; not every orientation was used in every patient. The shape, relative size, and obliquity of the simulations matched the measurements reasonably well in most cases. However, the measured focal heating was more diffuse than the simulations predicted. The orange lines indicate the 50% contours for the MRTI; yellow lines indicate 25% contours for the the difference between the predictions and measurements was constant. However, in 15/32 patients, the simulations began to overpredict the focal heating as the energy increased. These deviant sonications (red crosses) were excluded from the attenuation optimization. Using the optimized density/sound speed relationship resulted in better agreement for the deviant sonications (c). The lines in (d)-(f) connect data obtained in each patient. (S, YI: slope, y intercept of linear regression).

FIG. 8.
Comparing predicted and measured focal temperatures. (a) Measured and simulated focal heating as a function of the applied acoustic energy. For the sonications where the measurements and simulations deviated (red crosses in Fig. 7   Optimizing density/attenuation and density/sound speed relationships using results from 72 patients. The relationships found were similar for those for patients 1-32.

FIG. 13.
Impact of spatial averaging and aberration correction on focal heating. (a) Effects of spatial averaging on the simulated peak temperature rise at the focus using the phase corrections calculated by the manufacturer, the ideal phase corrections, and no correction. The models suggest that the peak temperature rise was 1.4-1.7 times the value measured with MRTI. (b) Plot of simulated temperature rise using ideal phase corrections and no correction versus those used during treatment (Tx). The models suggest that the correction used during the treatment improved the measured heating by a factor of 2, and that using the ideal correction predicted could yield a mean improvement of approximately 20% (dotted lines: linear regressions; solid lines: unity).

FIG. 14.
Skull-derived metrics and treatment efficiency. Maps of the outer skull surface color-coded based on different metrics derived from the CT scans are shown for 32 patients. The patients are sorted by the acoustic energy needed to reach a focal temperature of 55 °C, which covered an order of magnitude. A general trend towards increasing acoustic energies as the loss due to attenuation, skull thickness, and attenuation coefficient of the diploe increased and as the trabecular/cortical density ratio decreased is evident. Clear outliers for each metric are also evident.  Coefficients for Eqs. (3) and (6) that resulted in the best agreement between MRTI and the simulations.