Comparison between MR and CT imaging used to correct for skull

Scientific Reports volume 12, Article number: 13407 (2022) Cite this article

1764 Accesses

2 Citations

4 Altmetric

Metrics details

Transcranial focused ultrasound with the InSightec Exablate system uses thermal ablation for the treatment of movement and mood disorders and blood brain barrier disruption for tumor therapy. The system uses computed tomography (CT) images to calculate phase corrections that account for aberrations caused by the human skull. This work investigates whether magnetic resonance (MR) images can be used as an alternative to CT images to calculate phase corrections. Phase corrections were calculated using the gold standard hydrophone method and the standard of care InSightec ray tracing method. MR binary image mask, MR-simulated-CT (MRsimCT), and CT images of three ex vivo human skulls were supplied as inputs to the InSightec ray tracing method. The degassed ex vivo human skulls were sonicated with a 670 kHz hemispherical phased array transducer (InSightec Exablate 4000). 3D raster scans of the beam profiles were acquired using a hydrophone mounted on a 3-axis positioner system. Focal spots were evaluated using six metrics: pressure at the target, peak pressure, intensity at the target, peak intensity, positioning error, and focal spot volume. Targets at the geometric focus and 5 mm lateral to the geometric focus were investigated. There was no statistical difference between any of the metrics at either target using either MRsimCT or CT for phase aberration correction. As opposed to the MRsimCT, the use of CT images for aberration correction requires registration to the treatment day MR images; CT misregistration within a range of ± 2 degrees of rotation error along three dimensions was shown to reduce focal spot intensity by up to 9.4%. MRsimCT images used for phase aberration correction for the skull produce similar results as CT-based correction, while avoiding both CT to MR registration errors and unnecessary patient exposure to ionizing radiation.

Transcranial focused ultrasound (FUS) is an incisionless therapeutic modality that has human applications for thermal ablation1,2,3,4 as an alternative to open brain surgery, neuromodulation5,6,7 to probe the circuitry of the brain, and blood brain barrier opening8,9,10,11 to improve targeted drug delivery into the brain. In these applications, ultrasound is transmitted through the intact skull and focused deep into brain tissue to induce thermal or biomechanical effects. These effects can either be transient or irreversible depending on the parameters used during treatment.

Focusing ultrasound through the intact skull is challenging because the skull defocuses the ultrasound beam and shifts the focal spot away from the target. Correcting for the skull’s effects is highly complex due to the skull’s heterogeneity within and across patients, often varying in thickness, shape, size, and bone composition. Differences between skulls result in a broad range of transcranial ultrasound efficiencies across patients, such that the same applied power can lead to a four-fold range in temperature rise at the focal spot12.

Many methods have been proposed to account for the aberrating effects of the skull in order to refocus the focal spot. These methods include ray tracing13,14,15, finite difference time domain (FDTD)14,16,17,18,19,20,21,22,23,24,25,26, pseudo-spectral time domain (PSTD)27,28,29,30,31, and hybrid angular spectrum (HAS)12,32,33,34,35,36,37. Although these methods are different in theory and implementation, they have one similarity: they need to computationally model the skull to estimate and correct for its aberrating effects. For clinical ablation treatments, computed tomography (CT) images are used to create these computational skull models. However, the use of CT has several downsides. First, pre-operative CT images need to be registered to intra-operative magnetic resonance (MR) images, thus potentially leading to registration errors. Errors in registration could compromise phase corrections and lead to degradation of beamforming performance. Second, CT exposes patients to ionizing radiation, which is undesirable, for example, in pediatric populations. It is also undesirable for use with neuromodulation, as volunteers may be discouraged from participating in experiments needed to translate the technology for human use. Third, using CT images increases treatment workflow overhead, requiring separate imaging sessions to acquire all pre-operative images.

An alternative to CT is MR imaging, a safe, non-ionizing imaging modality. In recent years, developments in MR imaging of bone have been able to produce images with CT-like bone contrast38,39,40,41,42,43. These MR imaging techniques use ultrashort echo time (UTE), zero echo time (ZTE), or T1-weighted magnetization preparation rapid acquisition gradient echo (T1w MPRAGE) sequences to acquire bone signal before it decays rapidly due to T2* effects. T1w MPRAGE-derived skull models performed comparably to CT-derived skull models when simulating a single element transducer43, and UTE-derived skull models performed comparably to CT-derived skull models when simulating a hemispherical head transducer41,42. Furthermore, when using the standard of care InSightec ray tracing method for phase correction, UTE-derived binary images of the skull performed comparably with CT images44. In both Guo et al. and Miller et al., MR proton resonance shift thermometry (MR thermometry) was used to evaluate non-simulated experimental results. However, due to the method’s underlying resolution of 1.09 × 2.18 × 3 mm, it outputs data that averages the intensity of the focal spot and quantizes positioning errors into step sizes that are relatively large compared to a nominal focal spot size of 1.5 × 1.5 × 3 mm.

In this study, we directly compared the effectiveness of using CT versus MR-simulated-CT (MRsimCT) images for beamforming in an InSightec Exablate phased array transducer. We performed 3D hydrophone raster scans to measure the resulting focal spots and evaluated them using six metrics: pressure at the target, peak pressure, intensity at the target, peak intensity, positioning error, and focal spot volume. Additionally, we investigated the effects of CT misregistration on beamforming to emulate potential registration errors introduced during a clinical treatment. This provided a clearer comparison between the uses of CT and MRsimCT for treatment planning.

Three different image contrasts were generated using the acquired MR data: MRsimCT38, –log(short TE)39, and short TE–long TE38 (Fig. 1, Supplementary Figs. 1, 2). Qualitatively, all three contrasts preserved the overall skull shape, size, and thickness. Of the three options, MRsimCT achieved the most similar skull density ratio (SDR) as the CT, followed by –log(short TE), then short TE–long TE. SDR is a measure of skull homogeneity, with a value of 0 signifying low homogeneity and a value of 1 signifying high homogeneity. Changes in SDR compared to the CT suggest that bone composition information was lost.

CT contrast of the skull compared to various MR contrasts (skull A). Cortical and trabecular bone contrast is clearly depicted by CT and may be preserved by two of the three MR post-processing methods. The units for MRsimCT, –log(short TE), and short TE–long TE are not the same, thus different windowing and leveling were used. Scale bars and skull density ratios (SDRs) are shown at the bottom of each image. MRsimCT was the preferred choice based on bone contrast, minimal SDR change compared to CT, minimal background signal bias, and generalizable post-processing. Therefore, it was used for the remainder of this study.

Amongst the three MR contrasts of MRsimCT, –log(short TE), and short TE–long TE, MRsimCT was the preferred choice based on bone contrast, minimal SDR change compared to CT, and minimal background signal bias (Fig. 1, Supplementary Figs. 1–3). Background signal bias was due to RF transmit and receive nonuniformities, coil sensitivities, and/or coil shading. Although –log(short TE) also preserved some bone contrast, the remaining background bias would be a source of error when mapping from image values to acoustic properties. The residual background bias may be an artifact of the bias correction algorithm39 because it was observed in both ex vivo and in vivo images39. Furthermore, the bias correction algorithm required substantial parameter tuning to yield sufficient efficacy with the input images. Therefore, it may need to be tailored to each set of MR acquisition parameters. In contrast, the computation of short TE–long TE can be done rapidly, though it also exhibits residual background signal bias (Fig. 1, compare the anterior and posterior portions of the skull). Unlike the other two options, MRsimCT images exhibit minimal background signal bias likely because of its post-processing:

where MTEshort and MTElong are the magnitude images from the short and long echo times, STEshort and STElong are the tissue signals from the short and long echo times, and B is the spatially nonuniform background bias due to RF transmit and receive nonuniformities, coil sensitivities, and/or coil shading. The magnitude images are a result of the underlying tissue signal weighted by the background bias term. From the MRsimCT post-processing, the bias term can be factored out from the numerator and denominator and removed from the overall expression.

MRsimCT was used for the remainder of this study for calculating phase corrections using the InSightec algorithm. Because only three skulls were used in this study, a nominal linear relationship between CT and MRsimCT values was used to estimate CT Hounsfield units (HU) (Fig. 2). Linear regressions fit to the data39,41 may also be used, though they are suitable only when large amounts of data are available. Otherwise, the linear regression may overfit to the data and result in loss of generalizability.

Relationship between CT and MRsimCT values. MRsimCT values prior to HUbone scaling in (2) span a range from 0 to 1 and can be used as an estimate for bone fraction. For the CT parameters used in this study, pure bone was calculated to have a value of 2000 HU36. The black line depicts the nominal relationship between CT and MRsimCT values. Because only three ex vivo skulls were used in this study, using a linear regression to predict CT HU from MRsimCT values may overfit to the data. Therefore, the nominal linear relationship was used to predict CT HU instead.

3D hydrophone raster scans were acquired for five sets of phase corrections: no correction, InSightec MR binary image mask, InSightec MRsimCT, InSightec CT, and hydrophone. The experimental setup is shown in Fig. 3, and example normalized intensity images from skull A are shown in Fig. 4. Three-plane cross sections of the focal spots can be found in Supplementary Figs. 4–6. With no phase correction, the focal spot was defocused and displaced from the geometric focus. The phase correction methods were able to recover the focal spot at or near the geometric focus, with varying pressures and positions. Using the hydrophone method for comparison, the presence of the skull reduced transmitted pressure to 25 ± 3% compared to free field. The focal spots were evaluated using six metrics, which are shown in Fig. 5 and summarized in Table 1.

Experimental setup used in this study. Three ex vivo skulls were used for experimentation. Each of the three skulls was fixed to a head frame to ensure consistent positioning. The head frame was secured in an InSightec 670 kHz hemispherical phased array transducer. Data were acquired with a needle hydrophone and 3-axis positioner system. Illustrations were drawn by Sarah Hwang.

2D cross section raster scans of the refocused focal spot (skull A). The sagittal view is shown. These cross sections are a subset of the 3D volume scans acquired. A red x marks the position of the target, which was placed at the geometric focus. The focal spot maximum may be out of plane. Figure 5 shows quantitative metrics for all three skulls.

Beamforming performance of different image types with the target placed at the geometric focus. (a) Target and peak pressure (shown with darker and lighter colors, respectively) when operating the transducer at 20 W of electrical power. (b) Target and peak intensity normalized to the hydrophone method. (c) Focal spot positioning error. (d) Focal spot volume. (e) Dice similarity coefficient. Standard deviation bars are shown.

When targeting the geometric focus, the InSightec MR binary image mask, InSightec MRsimCT, and InSightec CT achieved comparable phase correction performance with respect to normalized target intensity, normalized peak intensity, positioning error, and focal spot volume (Fig. 5, Table 1). A one-way repeated measures ANOVA was performed, followed by a Tukey’s multiple comparisons test. No statistical difference was observed for any of the metrics (Supplementary Table 1). When targeting 5 mm to the left of the geometric focus, the InSightec MRsimCT and InSightec CT continued to perform similarly, whereas the InSightec MR binary image mask experienced degradation in performance. In particular, the target and peak intensities for InSightec MR binary image mask substantially decreased compared to the geometric focus case (Supplementary Table 2, Supplementary Fig. 8). A one-way repeated measures ANOVA and Tukey’s multiple comparisons test were performed, and there was again no statistical difference between image types for any of the metrics (Supplementary Table 3).

To investigate whether the nominal linear relationship between CT and MRsimCT introduced systematic beamforming errors, we additionally calculated phase corrections using a linear relationship fit to the data (Supplementary Fig. 9). Target pressures were computed using the post hoc analyses described in the methods. There was no statistical difference between the nominal linear relationship and the linear relationship fit to the data.

One benefit to using MR images is foregoing the need to perform image registration between CT and MR images. By registering pre-operative CT images with intra-operative MR images, misregistration errors can be introduced and reduce the effectiveness of estimated phase corrections. We performed error sensitivity analyses to estimate the influence of these registration errors. When errors were applied along one dimension, normalized target intensities decreased by approximately 10% within 4 mm of displacement (Fig. 6a) and within 4 degrees of rotation (Fig. 6b). 20% loss in normalized target intensity was observed within 6 mm of displacement and within 8 degrees of rotation. Phase correction performance was less sensitive to displacements along the superior axis, which may be because incidence angles were least affected by displacements along that axis. With rotation errors, performance was similarly sensitive along all three cardinal axes.

Effects of CT misregistration on normalized target intensity. (a) Displacement and (b) rotation errors were applied along one dimension only. Axes are colored according to their corresponding dimensions.

In practice, there exists an upper bound on the displacement and rotation errors due to misregistration. The upper bounds were approximated by performing manual registrations between pre-operative CT and intra-operative MR images, then evaluating the variability across multiple registrations. For two patients, manual registrations were performed five times each. Each registration was used as “ground truth” while the other four registrations were used to calculate the mean absolute registration errors. Across all patient registrations, the mean displacement errors were less than 0.4 mm and the mean rotation errors were less than 2 degrees. In keeping with the error analysis resolution of 1 mm and 1 degree for displacement and rotation, respectively, we assumed negligible displacement error and focused our analyses on combinations of rotation errors along three different axes. When rotation errors were applied in two dimensions, the maximum normalized target intensity loss within a ± 2 degree range was 8.8%. With rotation errors in three dimensions, the maximum normalized target intensity loss within a ± 2 degree range was 9.4%.

In this study, we compared the quality of several MR contrasts for bone imaging and directly compared the efficacy of several image types for correcting skull-induced phase aberrations. Although the MR binary image mask, MRsimCT, and CT images performed similarly when targeting the geometric focus, MR binary image mask performance substantially degraded when targeting 5 mm left of the geometric focus. MRsimCT shows great potential to be used as an alternative to CT. Qualitatively, the images exhibit CT-like bone contrast and preserve information on skull shape, size, thickness, and bone composition. Quantitatively, both sets of images resulted in similar target and peak intensities, positioning error, and focal spot volume. Investigation of CT misregistration showed that within the practical limits of error, the reduction in target intensity may be clinically relevant for ablation treatments.

The beamforming performances of InSightec MRsimCT and InSightec CT were comparable based on six quantitative metrics. The same standard of care InSightec phase correction algorithm was used on both types of images, and the results suggest that MRsimCT images may already be suitable to be integrated into the existing clinical workflow for treatment planning.

Miller et al. previously reported that an MR binary image mask performed as well as CT when used with the InSightec ray tracing method44. Even with 10 mm of steering away from the geometric focus, the measured MR thermometry peak temperature rise and focal spot position were highly comparable, with less than 2.3% and 5.8% difference, respectively. The results from the present study show larger differences between MR binary image mask and CT, with 22.4% difference in peak intensity and 49.6% difference in positioning error when steering 5 mm left of the geometric focus. The increase in position error difference observed in this study could be due to the higher spatial resolution of hydrophone raster scans (0.25 mm isotropic) compared to MR thermometry (typically 1.09 × 2.18 × 3 mm). The increase in peak intensity difference remains unclear, though it could be due to the use of different ex vivo skulls or InSightec software versions. Unfortunately, we were unable to determine whether the same skulls or software versions were used in both studies.

Post hoc analyses showed that there was no statistical difference in beamforming performance between the nominal linear relationship and the linear relationship fit to the data, however, there are a couple of caveats to using the latter. One caveat is the potential loss of generalizability when there is a small amount of available data for model training and testing. Another caveat is the linear regression’s dependence on the image acquisition parameters. Both CT and MR acquisition parameters need to be kept constant in order to maintain a stable relationship between their image values. As an example, a change in MR receive gain during acquisition could drastically affect the mapping from MR to CT values. Changes in parameters such as resolution, MR flip angle, or CT tube voltage45 would also have non-negligible effects. Challenges exist in maintaining a fixed set of image acquisition parameters across operators and across imaging sites; despite a CT imaging protocol established by InSightec, there still exists variability in CT acquisition parameters. Similarly, there will also be challenges in maintaining a fixed set of MR acquisition parameters across operators and across sites. Thus, for future studies, accounting for image acquisition parameters for both MR and CT will be of great importance.

Up until this point, MRsimCT and CT were compared based on their best-case performances. In this study, both sets of images were well registered to the hemispherical transducer due to the use of fiducial markers in each head holder. In contrast, during clinical treatments, pre-operative CT images need to be registered to intra-operative MR images without the use of fiducial markers. Furthermore, registration to the intra-operative MR images is additionally challenging because of B1 inhomogeneity artifacts created by the ultrasound hemispherical transducer and water surrounding the head. This process can introduce registration errors, which may reduce the effectiveness of estimated phase corrections. Therefore, applying displacement and rotation errors to the registered CT images (and not introducing any registration errors to the MRsimCT images) allowed for a fairer comparison between MRsimCT and CT.

Registration errors reduced the normalized target intensity by a non-negligible amount that may be clinically relevant for ablation treatments. Although there exists a large parameter space of misregistration errors that lead to large losses in intensity, in practice, the errors are likely bounded to within ± 2 degrees of rotation along each dimension. Within this range, the maximum normalized target intensity loss was 8.8% when errors were present in two dimensions and 9.4% when errors were in three dimensions. For ablation treatments, this loss in intensity may not drastically affect patients whose skulls permit efficient transmission of ultrasound. However, for patients with highly aberrating or attenuating skulls, the intensity loss could be the difference between successfully or unsuccessfully achieving the desired ablative temperature at the target. The limitation is often not transducer output power, but rather, increased patient discomfort when operating at high powers.

Although the MRsimCT versus CT results are encouraging, they should be considered with the limitations of this study. First, experiments were performed on only three skulls, a small sample size that may not be representative of the entire patient population. Investigation into additional skulls will need to be performed, ideally spanning a larger range of SDR values. Second, only two targets were used for phase correction. Although this may be sufficient for treatments requiring a small window of treatment, further investigation with electronic steering will be needed. Third, comparisons between the three types of images were performed using only the standard of care InSightec ray tracing algorithm. If a more accurate simulation method were used instead, there may be more pronounced differences in performance between the types of images. Although the conclusions presented in this study are applicable to the current standard of care, they may have to be revisited when considering new phase correction algorithms used in the future. These study limitations motivate important future work to acquire additional data and to characterize the effectiveness of MRsimCT versus CT for transcranial treatment planning.

Compared to the use of CT images for treatment planning, the use of MRsimCT images produced similar results. Using MR images as an alternative to CT images has great relevance for transcranial FUS ablation treatments. It would reduce beamforming errors due to misregistration and avoid patient exposure to ionizing radiation. These improvements are also highly relevant for applications in neuromodulation and blood brain barrier opening, which are currently being translated for human use.

Three ex vivo human skulls were used in this study. The protocol was approved by the Virginia State Anatomical Program and all research was conducted in accordance with the guidelines specified by the Virginia State Anatomical Program, the University of Virginia medical research protocols, and the University of Virginia Institutional Review Board for Health Sciences Research. Informed consent to use the donors’ skulls for scientific research were obtained from the donors and their next of kin.

The skulls were each fixed to a head frame prior to experimentation to ensure reproducible positioning between experiments (Fig. 3). Prior to imaging or other experimentation, the skulls were degassed overnight with an Abbess Instruments acrylic vacuum chamber. After degassing, the skulls were transferred to low density polyethylene bags in preparation for computed tomography (CT) and magnetic resonance (MR) imaging. The transfers were performed without exposing the skulls to air. A margin of 40 mm of degassed water between the skull and the bag material was used to avoid MR susceptibility artifacts.

Computed tomography (CT) scans for all three skulls were acquired using CT parameters approved by the InSightec patient screening imaging protocol (Table 2). The skull density ratios (SDR) are also reported. In general, SDR is a measure of bone composition homogeneity on a scale of 0 to 1, with higher SDR associated with higher homogeneity.

Ultrashort echo time (UTE) magnetic resonance (MR) images were acquired on a Siemens 3T Prisma MR scanner with a clinical 3T transmit-receive circularly polarized head coil (Siemens Healthcare GmbH). A 3D stack-of-spirals sequence46 was used with TE 50 μs and 2510 μs, TR 11 ms, flip angle 20 degrees, 1049 spiral interleaves, and matrix size 416 × 416 × 192 for a nominal resolution of 0.82 × 0.82 × 0.80 mm assuming no T2* decay. The UTE images were used to generate the MR-simulated-CT (MRsimCT) images38 using the following expression:

where MTEshort is the magnitude image from the short echo time, MTElong is the magnitude image from the long echo time, and HUbone is the maximum HU of bone. HUbone was calculated to be 2000 HU for the CT parameters used in this study using ρbone and mass attenuation coefficients for bone and water reported by the National Institute of Standards and Technology36,47 (Supplementary Fig. 10). We assumed an effective tube voltage of 60 kV for a 120 kVp spectrum because the spectrum data was not available to us. Prior to scaling by HUbone, the MRsimCT values spanned a range from 0 to 1 and can be used as an estimate for bone fraction. A value of 0 denotes water or a small proportion of bone in the voxel and a value of 1 denotes a large proportion of bone. Therefore, performing the HUbone scaling gives a nominal relationship with which to predict CT HU. Although a linear regression between CT HU and MRsimCT values could also have been used, it may overfit to the data because only three skulls were used in this study. MR binary image masks were generated using a threshold of 0 HU. Bone voxels were set to 1000 HU and all other voxels were set to − 1000 HU44. Creating the MR binary image masks as described by Miller et al.44 presented the opportunity to benchmark the results from this study against the results from that study. This was of interest because both studies were performed at the University of Virginia using the same InSightec equipment. As an alternative to the MR binary image mask, additional complexity may be added; a multi-layered skull approach would be an intermediate step between the MR binary image mask method and the MRsimCT method. As seen in Fig. 1, there is enough bone contrast to delineate between cortical and trabecular bone. Either nominal HU values or average HU values calculated from a cohort of patients could be assigned to the bone layers. In this study, to use results from Miller et al.44 as points of comparison, an MR binary image mask was used instead of a multi-layered image mask.

Two other MR contrasts were investigated in this study: –log(short TE)39 and short TE–long TE38. Prior to calculating –log(short TE), the spatially nonuniform background bias was removed from the magnitude images using a bias correction algorithm reported by Wiesinger et al.39 (for these data, K = 60 levels of resolution and a low pass filter three times larger than the smallest region of interest were used). For the short TE–long TE computation, background bias was not removed.

The skull and head frame were secured in a 670 kHz InSightec Exablate 4000 hemispherical phased array transducer and positioned in the clinical orientation; the anterior portion of the skull was closest to the anterior portion of the transducer. An Onda HNA-0400 needle hydrophone was mounted on an Onda 3-axis positioner system and used to perform 3D raster scans of the beam profile.

Prior to skull sonication, the position of the geometric focus in water was determined using a series of 2D raster scans. All transducer elements were turned on, and the focal spot position was localized by alternating between axial, coronal, and sagittal scan planes and centering the system on the position with maximum pressure. A three-plane cross section of the focal spot in water is shown in Supplementary Fig. 11.

The InSightec workstation was operated in research mode and was supplied phase corrections prior to each sonication. All elements were fired simultaneously while the hydrophone was rastered in a 5 × 5 × 5 mm field of view grid. This 3D grid was acquired with 21 individual 5 × 5 mm sagittal scans with in-plane and through-plane step sizes of 0.25 mm. For each skull, data for all 21 sagittal scans and for all phase correction schemes were acquired on the same day. The pulsing scheme was 2 ms on, 50 ms off using 20 W of electrical power. An Agilent DS07012B oscilloscope and the Soniq software from Onda were used to record the peak to peak voltages, which were then converted into pressure and intensity amplitudes.

To accurately model the experimental setup in the InSightec system, the position and orientation of the skull relative to the transducer was determined. Registration between the skull and head frame was not necessary, because their relative positioning and orientation were captured in the CT and MR images. The head frame was fabricated according to computer aided design (CAD) specifications, in which the position of the head frame relative to the transducer was also specified. Eight 2 mm diameter tantalum beads were epoxied onto the head frame at specified locations denoted in the CAD. These beads were installed prior to imaging and served as fiducial markers to register the images to the transducer (Supplementary Fig. 12). The centroid of each bead was used as its position in image space. A singular value decomposition-based least squares registration15,37,48 was performed to achieve point-wise registration between fiducial marker positions in image and transducer space. The transformation matrix that registers the two coordinate systems is:

where Ptransducer is a 3 × 8 matrix of fiducial marker positions in transducer space, Pimage is a 3 × 8 matrix of fiducial marker positions in image space, R is a 3 × 3 rotation matrix, and Tr is a 3 × 8 translation matrix. The rotation matrix R is calculated as

where U and V are the unitary matrices from the singular value decomposition of the matrix H

and where pimage,i is a 3 × 1 position vector of the ith fiducial marker in image space, ptransducer,i is a 3 × 1 position vector of the ith fiducial marker in transducer space, and T is the transpose operator. With values for Pimage, Ptranducer, and R, the translation matrix Tr can then be calculated using (3). The resulting transformation matrix was used to resample the images into transducer space to create a set of registered images. Axial slices with 0.5 × 0.5 × 1 mm resolution were used for the MR binary image mask, MRsimCT, and CT images.

The geometric focus and 5 mm left of the geometric focus were selected as targets, and four different methods were used to estimate phase corrections. These methods were (1) InSightec ray tracing with MR binary image mask inputs, (2) InSightec ray tracing with MRsimCT image inputs, (3) InSightec ray tracing with CT image inputs, and (4) hydrophone.

The InSightec ray tracing method requires computational models of the skulls to compute phase corrections. These skull models were generated using either MR binary image mask, MRsimCT, or CT images of the skulls, which were first registered to the transducer as described above.

Although the InSightec ray tracing method remains proprietary, we were able to calculate phase corrections by using the InSightec workstation in clinical mode. After uploading the registered MR binary image mask, MRsimCT, or CT images, we followed the clinical workflow for treatment planning (Supplementary Fig. 13). Because the skull was already registered to the transducer in these images, no further registration was necessary. Under the hood, the InSightec workstation maps image voxel values to acoustic properties, which are then used to calculate phase corrections. The phase corrections could be found in the system logs.

While the InSightec method requires models of the skull to compute phase corrections, the hydrophone method directly accounts for skull aberrations using experimental data. Individual elements were fired sequentially, and their time series signals were measured with the hydrophone positioned at the geometric focus. The pulse duration used was 3 ms at 375 W electrical power, the sampling frequency was 15 MHz, and the sampling duration was 7 ms. Data were acquired for all 1024 transducer elements using a TiePie Handyscope HS6 Diff and a Panametrics model 5676 preamplifier.

Six quantitative metrics were used to evaluate phase correction performance. Focal spot pressures were indexed at the target position and at the position of peak pressure. Intensities were calculated and normalized to the results from the hydrophone phase correction method, the gold standard for recovering maximum intensity. This normalization removed the influence of skull-dependent attenuation and allowed phase correction performance to be compared across skulls. Positioning errors were calculated as the difference between the target position and the focal spot position, defined to be at the voxel of peak pressure. Focal spot volume was measured at pressure full width half maximum.

Post hoc analyses of the (i) MRsimCT to CT relationship and (ii) CT misregistration were performed. To perform these analyses, phase corrections were applied to single element time series data previously acquired with the hydrophone method, detailed in the “Phase correction methods” section. By the superposition principle, the net focal spot pressure is equivalent to the sum of pressures contributed by the individual elements. Because the single element data were only acquired at the target position, only the target pressure could be computed for these analyses. The target pressure resulting from any set of phase corrections can be determined by the following expression:

where P is the pressure averaged over 100 cycles, si is the time series data for element i, and \({\phi}\)i is the phase correction for element i calculated with one of the methods. Analyses by this method removed the need to reacquire raster scan data for each set of new phase corrections and substantially reduced the number of data acquisitions needed.

To investigate the effects of CT misregistration, displacements and rotations were applied to the registered CT images prior to computing phase corrections. Misregistration errors can cover a large parameter space due to six degrees of freedom: three directions of displacement and three axes of rotation. Therefore, to reduce the complexity of analyses and data visualization, we first evaluated phase correction performance with errors applied in one dimension. Displacement errors were applied in 1 mm increments along the left–right, posterior-anterior, or superior-inferior axes. Rotation errors were applied along the same axes, with the origin centered at the focal spot target. Increments of 1 degree were applied up until ± 10 degrees, after which increments of 5 degrees were applied. Phase corrections were recalculated for each displaced or rotated skull model. Post hoc analyses were then performed to evaluate the effects of misregistration on the resulting target pressure.

Elias, W. J. et al. A pilot study of focused ultrasound thalamotomy for essential tremor. N. Engl. J. Med. 369, 640–648 (2013).

Article CAS PubMed Google Scholar

Coluccia, D. et al. First noninvasive thermal ablation of a brain tumor with MR-guided focused ultrasound. J. Ther. Ultrasound 2, 17 (2014).

Article PubMed PubMed Central Google Scholar

Magara, A. et al. First experience with MR-guided focused ultrasound in the treatment of Parkinson’s disease. J. Ther. Ultrasound 2, 11 (2014).

Article PubMed PubMed Central Google Scholar

Martínez-Fernández, R. et al. Focused ultrasound subthalamotomy in patients with asymmetric Parkinson’s disease: A pilot study. Lancet Neurol. 17, 54 (2018).

Article PubMed Google Scholar

Legon, W. et al. Transcranial focused ultrasound modulates the activity of primary somatosensory cortex in humans. Nat. Neurosci. 17, 322–329 (2014).

Article CAS PubMed Google Scholar

Lee, W. et al. Image-guided transcranial focused ultrasound stimulates human primary somatosensory cortex. Sci. Rep. 5, 8743 (2015).

Article PubMed PubMed Central CAS Google Scholar

Lee, W. et al. Transcranial focused ultrasound stimulation of human primary visual cortex. Sci. Rep. 6, 34026 (2016).

Article ADS CAS PubMed PubMed Central Google Scholar

Patrick, J. T. et al. Ultrasound and the blood-brain barrier. In Consensus on Hyperthermia for the 1990s (eds Bicher, H. I. et al.) (Springer, 1990).

Google Scholar

Hynynen, K. & Mcdannold, N. Imaging-guided focal opening of the blood-brain barrier in rabbits. Radiology 220, 640–646 (2001).

Article CAS PubMed Google Scholar

Mcdannold, N., Arvanitis, C. D., Vykhodtseva, N. & Livingstone, M. S. Temporary disruption of the blood–brain barrier by use of ultrasound and microbubbles: Safety and efficacy evaluation in Rhesus Macaques. Am. Assoc. Cancer Res. 72, 3652–3664 (2012).

Article CAS Google Scholar

Downs, M. E., Buch, A., Karakatsani, M. E., Konofagou, E. E. & Ferrera, V. P. Blood-brain barrier opening in behaving non-human primates via focused ultrasound with systemically administered microbubbles. Sci. Rep. 5, 15076 (2015).

Article ADS CAS PubMed PubMed Central Google Scholar

Vyas, U. et al. Predicting variation in subject thermal response during transcranial magnetic resonance guided focused ultrasound surgery: Comparison in seventeen subject datasets. Med. Phys. 43, 5170–5180 (2016).

Article PubMed PubMed Central Google Scholar

Kyriakou, A. et al. A review of numerical and experimental compensation techniques for skull-induced phase aberrations in transcranial focused ultrasound. Int. J. Hyperth. 30, 36–46 (2014).

Article Google Scholar

Kyriakou, A., Neufeld, E., Werner, B., Székely, G. & Kuster, N. Full-wave acoustic and thermal modeling of transcranial ultrasound propagation and investigation of skull-induced aberration correction techniques: a feasibility study. J. Ther. Ultrasound 3, 11 (2015).

Article PubMed PubMed Central Google Scholar

Jin, C., Moore, D., Snell, J. & Paeng, D. An open-source phase correction toolkit for transcranial focused ultrasound. BMC Biomed. Eng. 2, 1–11 (2020).

Article Google Scholar

Connor, C., Clement, G. & Hynynen, K. A unified model for the speed of sound in cranial bone based on genetic algorithm optimization. Phys. Med. Biol. 47, 3925–3944 (2002).

Article PubMed Google Scholar

Aubry, J. F., Tanter, M., Pernot, M., Thomas, J. L. & Fink, M. Experimental demonstration of noninvasive transskull adaptive focusing based on prior computed tomography scans. J. Acoust. Soc. Am. 113, 84–93 (2003).

Article ADS CAS PubMed Google Scholar

Jones, R. M. & Hynynen, K. Comparison of analytical and numerical approaches for CT-based aberration correction in transcranial passive acoustic imaging. Phys. Med. Biol. 61, 23–36 (2016).

Article PubMed Google Scholar

Marsac, L. et al. Ex vivo optimisation of a heterogeneous speed of sound model of the human skull for non-invasive transcranial focused ultrasound at 1 MHz. Int. J. Hyperth. 33, 6350 (2017).

Article Google Scholar

Marquet, F. et al. Non-invasive transcranial ultrasound therapy based on a 3D CT scan: Protocol validation and in vitro results. Phys. Med. Biol. 54, 2597–2613 (2009).

Article CAS PubMed Google Scholar

Deffieux, T. & Konofagou, E. Numerical study of a simple transcranial focused ultrasound system applied to blood-brain barrier opening. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 57, 2637–2653 (2010).

Article PubMed PubMed Central Google Scholar

Pichardo, S. & Hynynen, K. Multi-frequency characterization of speed of sound for longitudinal transmission on freshly excised human skulls. Phys. Med. Biol. 1215, 282–286 (2011).

Google Scholar

Leduc, N., Okita, K., Sugiyama, K., Takagi, S. & Matsumoto, Y. Focus control in HIFU therapy assisted by time-reversal simulation with an iterative procedure for hot spot elimination. J. Biomech. Sci. Eng. 7, 43–56 (2012).

Article Google Scholar

Bouchoux, G., Shivashankar, R., Abruzzo, T. A. & Holland, C. K. In silico study of low-frequency transcranial ultrasound fields in acute ischemic stroke patients. Ultrasound Med. Biol. 40, 1154–1166 (2014).

Article PubMed PubMed Central Google Scholar

Pulkkinen, A., Werner, B., Martin, E. & Hynynen, K. Numerical simulations of clinical focused ultrasound functional neurosurgery. Phys. Med. Biol. 59, 1679–1700 (2014).

Article PubMed PubMed Central Google Scholar

Jones, R. M. & Reilly, M. A. O. Experimental demonstration of passive acoustic imaging in the human skull cavity using CT-based aberration corrections. Med. Phys. 42, 4385 (2015).

Article PubMed PubMed Central Google Scholar

Top, C. B., White, P. J. & McDannold, N. J. Nonthermal ablation of deep brain targets: A simulation study on a large animal model. Med. Phys. 43, 870–882 (2016).

Article PubMed PubMed Central Google Scholar

Robertson, J. L. B., Cox, B. T., Jaros, J. & Treeby, B. E. Accurate simulation of transcranial ultrasound propagation for ultrasonic neuromodulation and stimulation. J. Acoust. Soc. Am. 141, 1726–1738 (2017).

Article ADS CAS PubMed Google Scholar

Robertson, J., Martin, E., Cox, B. & Treeby, B. E. Sensitivity of simulated transcranial ultrasound fields to acoustic medium property maps. Phys. Med. Biol. 62, 2559 (2017).

Article PubMed Google Scholar

Mcdannold, N., White, P. J. & Cosgrove, R. Elementwise approach for simulating transcranial MRI-guided focused ultrasound thermal ablation. Phys. Rev. Res. 1, 033205 (2019).

Article CAS PubMed PubMed Central Google Scholar

Maimbourg, G. et al. Computationally efficient transcranial ultrasonic focusing: Taking advantage of the high correlation length of the human skull. IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 67, 1993 (2020).

Article Google Scholar

Vyas, U. & Christensen, D. Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method. IEEE Trans. Ultrason. Ferroelectr. Freq. Control 59, 1093–1100 (2012).

Article PubMed Google Scholar

Almquist, S., de Bever, J., Merrill, R., Parker, D. & Christensen, D. A full-wave phase aberration correction method for transcranial high-intensity focused ultrasound brain therapies. Annu. Int. Conf. IEEE Eng. Med. Biol. Soc. 2014, 310–313 (2014).

PubMed Google Scholar

Almquist, S., Parker, D. L. & Christensen, D. A. Simulation of hemispherical transducers for transcranial HIFU treatments using the hybrid angular spectrum approach. J. Ther. Ultrasound 3, P9 (2015).

Article PubMed Central Google Scholar

Almquist, S., Parker, D. L. & Christensen, D. A. Rapid full-wave phase aberration correction method for transcranial high-intensity focused ultrasound therapies. J. Ther. Ultrasound 4, 30 (2016).

Article PubMed PubMed Central Google Scholar

Leung, S. A., Webb, T. D., Bitton, R. R., Ghanouni, P. & Pauly, K. B. A rapid beam simulation framework for transcranial focused ultrasound. Sci. Rep. https://doi.org/10.1038/s41598-019-43775-6 (2019).

Article PubMed PubMed Central Google Scholar

Leung, S. A. et al. Transcranial focused ultrasound phase correction using the hybrid angular spectrum method. Sci. Rep. https://doi.org/10.1038/s41598-021-85535-5 (2021).

Article PubMed PubMed Central Google Scholar

Johnson, E. M., Vyas, U., Ghanouni, P., Pauly, K. B. & Pauly, J. M. Improved cortical bone specificity in UTE MR Imaging. Magn. Reson. Med. https://doi.org/10.1002/mrm.26160 (2016).

Article PubMed PubMed Central Google Scholar

Wiesinger, F. et al. Zero TE MR bone imaging in the head. Magn. Reson. Med. 114, 107–114 (2016).

Article Google Scholar

Wiesinger, F., Bylund, M., Yang, J. & Kaushik, S. Magnetic resonance in medicine zero TE-based pseudo-CT image conversion in the head and its application in PET/MR attenuation correction and MR-guided radiation therapy planning. Magn. Reson. Med. 80, 1440–1451. https://doi.org/10.1002/mrm.27134 (2018).

Article PubMed Google Scholar

Guo, S. et al. Feasibility of ultrashort echo time images using full-wave acoustic and thermal modeling for transcranial MRI-guided focused ultrasound (tcMRgFUS) planning. Phys. Med. Biol. https://doi.org/10.1088/1361-6560/ab12f7 (2019).

Article PubMed PubMed Central Google Scholar

Su, P., Guo, S., Roys, S., Maier, F. & Bhat, H. Transcranial MR imaging-guided focused ultrasound interventions using deep learning synthesized CT. Am. J. Neuroradiol. 41, 1–8 (2020).

Article Google Scholar

Koh, H., Park, T. Y., Chung, Y. A., Lee, J. & Kim, H. Acoustic simulation for transcranial focused ultrasound using GAN-based synthetic CT. IEEE J. Biomed. Health Inform. 26, 161–171 (2022).

Article PubMed Google Scholar

Miller, G. W., Eames, M., Snell, J. & Aubry, J. F. Ultrashort echo-time MRI versus CT for skull aberration correction in MR-guided transcranial focused ultrasound: In vitro comparison on human calvaria. Med. Phys. 42, 2223–2233 (2015).

Article PubMed Google Scholar

Webb, T. D. et al. 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, 1–20 (2018).

Article Google Scholar

Mugler, J. P. et al. Breath-Hold UTE Lung Imaging Using a Stack-of-Spirals Acquisition. (International Society for Magnetic Resonance Imaging, 2015).

Hubbell, J. H. & Seltzer, S. M. Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients. (National Institute of Standard Technology, 2004).

Arun, K., Huang, T. & Blostein, S. Least-squares fitting of two 3D point sets. IEEE Trans. Pattern Anal. Mach. Intell. 9, 698–700 (1987).

Article CAS PubMed Google Scholar

Download references

The authors would like to thank the Focused Ultrasound Foundation for their support, Nathaniel Kelm for his assistance with the InSightec hardware, and Sarah Hwang for her illustrations used in this article.

Funding was provided by the National Institutes of Health (R01 MH111825 and R01 EB028773).

Department of Bioengineering, Stanford University, Stanford, CA, USA

Steven A. Leung & Kim Butts Pauly

Focused Ultrasound Foundation, Charlottesville, VA, USA

David Moore & John Snell

Department of Biomedical Engineering, University of Virginia, Charlottesville, VA, USA

Yekaterina Gilbo, Craig H. Meyer & G. Wilson Miller

Department of Neurological Surgery, University of Virginia, Charlottesville, VA, USA

John Snell

Department of Electrical Engineering, Stanford University, Stanford, CA, USA

Taylor D. Webb & Kim Butts Pauly

Department of Radiology and Medical Imaging, University of Virginia, Charlottesville, VA, USA

Craig H. Meyer & G. Wilson Miller

Department of Radiology, Stanford University, Stanford, CA, USA

Pejman Ghanouni & Kim Butts Pauly

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

S.A.L. designed the experiments, analyzed the data, and prepared the manuscript. S.A.L. and D.M. acquired the ultrasound data. S.A.L., D.M., Y.G., C.H.M., and G.W.M. acquired the MR data. J.S. provided assistance with the data analysis. T.D.W., P.G., and K.B.P. provided many valuable discussions. K.B.P. supervised the research. All authors reviewed the manuscript.

Correspondence to Steven A. Leung.

The authors declare no competing interests.

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Leung, S.A., Moore, D., Gilbo, Y. et al. Comparison between MR and CT imaging used to correct for skull-induced phase aberrations during transcranial focused ultrasound. Sci Rep 12, 13407 (2022). https://doi.org/10.1038/s41598-022-17319-4

Download citation

Received: 15 November 2021

Accepted: 25 July 2022

Published: 04 August 2022

DOI: https://doi.org/10.1038/s41598-022-17319-4

Anyone you share the following link with will be able to read this content:

Sorry, a shareable link is not currently available for this article.

Provided by the Springer Nature SharedIt content-sharing initiative

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.