Comparison of transcranial MR-guided focused ultrasound phase correction simulations to clinically measured MRTI

Essential tremor and Parkinson's disease are neurological disorders causing involuntary shaking that affects an estimated 8 million people in the United States. Transcranial Magnetic Resonance-guided Focused Ultrasound (tcMRgFUS) is a noninvasive alternative method of treating tremors to the invasive deep brain stimulation standard treatment. Phase aberration correction is an important component of tcMRgFUS that best focuses the ultrasound beams into the brain. Simulations can be used to characterize the application of phase correction before starting the clinical treatment. They can also be used to compare different methods of phase aberration correction. In this study, three different methods of phase correction and the Magnetic Resonance Thermal Imaging (MRTI) from clinical treatments are created and compared. The simulation was performed by taking the Computed Tomography (CT) scan of a patient's head from clinical treatments and transforming it into MR treatment space. The CT is segmented and run through an acoustic and thermal simulation to simulate temperature rise from tcMRgFUS treatments. We found that the hottest voxel temperatures from clinical phase correction and MRTI had no significance differences for 16 out of 20 simulated patients. The other two methods of phase correction were time reversal and no phase correction. Time reversal on average produced higher temperatures and was significantly different than the other methods. Utilizing a method with no correction was significantly worse than the other methods of phase correction. Based on these findings, we will continue to simulate all clinical treatment patients.

COMPARISON OF TRANSCRANIAL MR-GUIDED FOCUSED ULTRASOUND PHASE CORRECTION SIMULATIONS TO CLINICALLY MEASURED MRTI by Emma N. Slominski A Senior Honors Thesis Submitted to the Faculty of The University of Utah In Partial Fulfillment of the Requirements for the Honors Degree in Bachelor of Science In Biomedical Engineering Approved: Dennis L. Parker, PhD Thesis Faculty Supervisor David W. Grainger, PhD Chair, Department of Biomedical Engineering Kelly W. Broadhead, PhD Honors Faculty Advisor Sylvia D. Torti, PhD Dean, Honors College May 2021 Copyright © 2021 All Rights Reserved ABSTRACT Essential tremor and Parkinson’s disease are neurological disorders causing involuntary shaking that affects an estimated 8 million people in the United States. Transcranial Magnetic Resonance-guided Focused Ultrasound (tcMRgFUS) is a noninvasive alternative method of treating tremors to the invasive deep brain stimulation standard treatment. Phase aberration correction is an important component of tcMRgFUS that best focuses the ultrasound beams into the brain. Simulations can be used to characterize the application of phase correction before starting the clinical treatment. They can also be used to compare different methods of phase aberration correction. In this study, three different methods of phase correction and the Magnetic Resonance Thermal Imaging (MRTI) from clinical treatments are created and compared. The simulation was performed by taking the Computed Tomography (CT) scan of a patient’s head from clinical treatments and transforming it into MR treatment space. The CT is segmented and run through an acoustic and thermal simulation to simulate temperature rise from tcMRgFUS treatments. We found that the hottest voxel temperatures from clinical phase correction and MRTI had no significance differences for 16 out of 20 simulated patients. The other two methods of phase correction were time reversal and no phase correction. Time reversal on average produced higher temperatures and was significantly different than the other methods. Utilizing a method with no correction was significantly worse than the other methods of phase correction. Based on these findings, we will continue to simulate all clinical treatment patients. ii 3 TABLE OF CONTENTS ABSTRACT ii INTRODUCTION 1 BACKGROUND 5 METHODS 11 RESULTS 16 DISCUSSION 23 REFERENCES 29 iii 1 INTRODUCTION Tremors from essential tremor and Parkinson’s disease cause involuntary, rhythmic muscle contractions [1]. These contractions primarily lead to shaking in the hands, but can also affect the arms, legs, torso, head, or even the vocal cords [1]. Essential tremor and Parkinson’s tremor both originate in the same location and can be treated in the same way [2]. However, essential tremor is eight times more common than Parkinson’s disease and affects 4% of adults over 40 years old [2], [3]. These neurological disorders progress with age and makes simple tasks like eating and writing difficult to impossible. Current treatments are limited to medication, therapy, and surgery. The current standard treatment for essential tremor uses invasive neurosurgery like deep brain stimulation [4]. However, invasive surgery is not a viable option for many elderly people due to health and safety risks A new alternative Federal Drug Administration (FDA) approved method of surgery is now available called transcranial magnetic resonance-guided focused ultrasound (tcMRgFUS). Magnetic resonance-guided focused ultrasound (MRgFUS) is a versatile and completely noninvasive technique that can perform targeted tissue thermal ablation to treat soft-tissue disorders such as uterine fibroids and breast tumors [5]. Transcranial MRgFUS is used to perform targeted tissue ablation in the brain. A common focus of MRgFUS simulation research addresses complications of transmitting acoustic energy through the hard tissue (i.e., skull bone) to access specific brain soft tissue regions needing ablation. Directed into the brain, MRgFUS can address soft tissue diseases such as brain tumors, essential tremor (ET), and Parkinson’s tremor [6]. This new method provides a totally noninvasive treatment option for patients not eligible for invasive neurosurgery. 2 Transcranial MRgFUS for essential tremor and Parkinson’s tremor has a set of criteria required in order to be eligible for the treatment. As the treatment involves magnetic resonance imaging (MRI), patients must not have metallic implants and cannot be claustrophobic. Another limiting factor is the skull density ratio (SDR) of a patient which is considered an indicator as to whether a tcMRgFUS treatment will be successful. The SDR assesses the skulls transparency to ultrasound, and low values are believed to be associated with lower treatment effectiveness [7]. Patients with an SDR<40 are generally said to be poor candidates for the procedure as more energy is required to ablate the tissue where tremors originate. More energy can also cause more damage to the skull and surrounding tissues. Better ways to predict treatment effectiveness outside of the SDR would be useful in determining if a patient is eligible for the procedure. Transcranial MRgFUS procedures require ultrasound energy from a large hemispherical ultrasound transducer external to the patient. A highly focused ultrasound beam through the intact skull bone produces an energy zone intended to ablate selected tissue inside the brain. Figure 1 depicts a view of focused ultrasound beams entering the head from the transducer. Due to individual variations in skull structure (e.g., thickness, density), the ultrasound energy transmission properties through the skull change from patient to patient [8]. Spatial bone heterogeneities and variations in skull thickness and composition of a patient’s skull affect the ultrasound beam path, resulting in focal spot position shifts, focal spot blurring, or the skull itself heating [9]. This departure from the ideal beam path is known as phase aberration, One way these issues have been avoided is by excluding patients with low SDRs to prevent performing treatments that are unlikely to be effective. Another way to address these issues is through phase aberration correction 3 algorithms that are applied to tcMRgFUS treatments to create a highly focused ablated lesion in each patient [10]. Phase correction can be calculated in many different ways, and acoustics simulations can be used to implement phase corrected tcMRgFUS in patient models [10], [11]. Figure 1: Depiction of ultrasound beams entering the skull. Three red simulated ultrasound beam waves move from the transducer elements to the focal spot in the brain. These beams are phase corrected and in phase at the focus after passing through the skull. Being in phase at the focus allows for the energy from each ultrasound beam to sum to a higher total energy at the focus. Skull is a sagittal side view. There are multiple different ways to perform phase aberration correction, however only one method of phase correction is approved for use in clinical treatments [4]. The method used in clinical treatments is a simple ray tracing method. Ray tracing functions by tracing each transducer element to the focal location, and the skull profiles along each ray are used to calculate the phase correction [12]. The absence of any phase correction is called no phase correction where ultrasound elements have a phase of zero. Another 4 method used in other MRgFUS applications is time reversal phase correction which backward propagates from a focal point source back to the transducer to create the phases based on a model [13]. This method could potentially be applied in transcranial applications. Time reversal, no phase correction, and the clinical ray tracing phase correction will all be analyzed in this research. We aimed to compare the different methods of phase correction in tcMRgFUS in thermal simulation modeling and to also compare it to clinical Magnetic Resonance Thermal Imaging (MRTI) that reports tissue heating data in patients for validation. By applying the simulation to match the MRTI scans of 20 patients, we are seeking the improved ability to be able to adjust the acoustic properties of the skull to find the model that best matches the heating observed in all patients [8]. To do so, we created tcMRgFUS simulations based on clinical data and implemented time reversal, no phase correction, and the clinical phase correction in these simulations to address heating. Different patients with their respective tcMRgFUS treatments were used to demonstrate how skull heterogeneities affect heating and phase correction. Producing and assessing thermal simulations was found to be useful to understand how effective a tcMRgFUS treatment might be in each specific patient before performing the treatment [14]. tcMRgFUS treatment effectiveness may be improved by simulating how the treatment would impact a patient model, and thereby knowing which phase correction method works best for each patient to ablate their soft tissue lesion. 5 BACKGROUND Acoustic Physics in Skull Tissues Acoustic properties affect how ultrasonic energy moves through bone and into tissue. As ultrasound propagates, the series of media in the head have different material properties and shapes, that attenuate, refract and distort the sound waves [8]. These heterogeneities defocus the energy and can create undesired heating. The interaction of the ultrasound with different layers of materials (e.g., air, bone, soft tissue) can be simulated to assess different contributions and produce corrections. One convenient method to produce the necessary properties and parameters for simulating ultrasonic energy traversing the skull for soft tissue ablations utilizes Computed Tomography (CT). CT images are reconstructed using multiple 2-D x-ray projection images to create a 3-D view of the body. CT images depict more attenuating tissues like bone as white and soft tissues like the brain as gray. The gradient of CT can be used to assign different properties to different tissues. CT and MRI scans are often taken before treatments to assist with diagnosis, as part of clinical studies, and help with treatment planning [15]. The CT scan taken before transcranial MRgFUS is used in the simulation for treatment planning to assess the SDR. One method to simulate the treatment is by segmenting the CT using thresholding. Thresholding involves assigning a number to a range of CT Hounsfield units. A Hounsfield unit (HU) is a linear transformation of the measured attenuation coefficient of the tissue [16]. By thresholding, similar tissues in the CT are grouped by their Hounsfield units and can be assigned acoustic properties. 6 The Hounsfield units in CT are used to separate the CT image volume into the different tissues of interest. Once separated, acoustic values are assigned to each type of tissue. The CT segmentation is a numeric model that splits different tissue types into groups. The segmented layers for cortical bone and diploe, i.e. cancellous bone, provide the biggest challenge for selecting values as there is a wide range of attenuation values for bone. Using incorrect attenuation values for bone could drastically change phase aberration correction calculations [17]. Phase aberration correction can be used to counter phase changes resulting from hard tissues. Phase aberration correction changes the phase angle so that ultrasound beams can focus within tissues based on the tissues they will be passing through [11]. Phase aberration correction changes the angle of each of the 1024 individual element beams in a transcranial focused ultrasound transducer. A focused ultrasound transducer creates ultrasound waves that shoot into a focused point. Hard tissues can change the direction of the ultrasound waves so that a focused point is no longer achieved [9]. By changing the phases, each beam’s phase error is corrected so that the maximum total intensity can be achieved in the focal spot in the brain. Acoustic properties dictate the beam distortion [18], phase correction causes the energy in the beams to combine constructively at the desired focus [19]. Phase correction is important because not accounting for heterogeneities in the head can drastically reduce the effectiveness of a treatment. Phase correction is created based on acoustic properties in numerical simulations to mimic the best possible phases for producing efficient ultrasound energy transfer through the skull and to a specific desired tissue focus for treatment. Phase correction is necessary because inhomogeneous tissues such as the skull bone can cause phase aberrations and shift the focal location [10]. Further, 7 the middle skull layer (diploe) can have structure that scatters the ultrasound energy and decreases the transmitted intensity. Patient selection for essential tremor treatments is further performed using a measurement called the skull-density-ratio (SDR), which utilizes the ratio of the cortical table and the diploe as measured in CT images. Patients with an SDR < 0.45 are generally excluded from treatments [7], [20] as skulls with low SDR can absorb most of the energy, injuring the patient skull while having an unsuccessful treatment [21]. As there is no other method used to predict heating before a treatment occurs, SDR has been used as the primary indicator as to whether a treatment will be successful or not. Treatments are reliant on SDR and the current clinical method of phase correction. It can be hypothesized that, with improved phase aberration correction methods, more of the ultrasound power transmitted from the transducer external to the patient can reach the intended focus at the tissue site, and more patients can be treated more effectively and safely. Phase correction is therefore implemented into tcMRgFUS treatments by adjusting the phases of a multi-element transcranial focused ultrasound transducer. tcMRgFUS Treatment Procedure The Insightec (Tirat Carmel, Israel) Exablate Neuro system is currently the only clinical tcMRgFUS system available for neurological applications [22]. Phase aberration correcting in the Exablate system is performed using a simple ray tracing method. Ray tracing finds all rays within each element and averages the phase changes over these rays [23]. With large-aperture transducers, such as the 30-cm hemispherical transducer in the Exablate Neuro system, phase aberration is more prominent than in other treatments using 8 smaller transducers due to the greater differences in acoustic properties experienced by the ultrasound beams emitted from the many elements in the larger transducer. TcMRgFUS treatments are performed on patients inside a 3-T Magnetic Resonance imaging (MRI) scanner with the Exablate device transducer. The transducer used in this research contains 1,024 elements in order to evenly distribute the ultrasound over as much of the head as possible. This is to assist in preventing the skull itself from becoming heated and damaged. Another way of preventing damage to the patient is by shaving their head to minimize energy absorption. Additionally, water (15-20°C) is filled into the volume between the bare scalp and the transducer to cool the scalp [4]. Once the patient is in position in the transducer, low power sonications are applied anywhere from 10-20 seconds until a peak temperature of 40-42°C is reached at the focus. These low powered sonication applications are then replaced with ones at higher power after it is determined through MRI and MR thermometry as well as patient feedback that the target area is being ablated. The high power sonication applications increase until the focal spot reaches 54°C or above. This temperature assures desired target tissue thermal ablation and lesion formation (necrosis), and the treatment is concluded [4]. MRI and sonication data from treatments like this are what will be used in this study. In the study, CT scans from patients are used to provide inputs for patient-specific simulations. Published Methods to Address De-Focusing One current approach used in MRgFUS simulations is the Hybrid Angular Spectrum (HAS) algorithm [22] – [24]. HAS has been modified for use in transcranial applications by replicating the hemispherical shape of the transducer [25]. Therefore, 9 simulations can be performed for brain treatments using this software. The software uses acoustic properties based on a segmented image from CT to create a pressure pattern (W/m3) at the same resolution as the segmented CT used. The pressure pattern is a result of simulated ultrasound beams propagating through the segmented image. The ultrasound beams can propagate using different methods of phase correction that include time reversal, ray tracing, implementing pre-calculated beams, and using no phase correction [22]. The Insightec system for current treatments uses a method of ray tracing for its phase aberration correction in patients. HAS is a crucial software tool to create acoustic simulations. Another method of creating acoustic simulations is through the use of k-Wave, an open source MATLAB toolbox that was utilized in the McDannold paper [21]. The McDannold paper uses many methods similar to what was used in this study; however the focus of their research was to create thermal simulations of patient data to predict where damage in the skull might occur. This differs from the purpose of this paper as we are looking to predict the lesion in the brain, and not the heating in the skull. However, the heating in the skull will still be observable. The skull heating location was compared and analyzed for how well they matched to clinical treatments. For the study in this paper, the focal point heating of the focused ultrasound is what will be primarily compared. Leung et al. [26] used the HAS method for acoustic simulations and also focused on looking at focal spot heating and position. The output from HAS was not specified, but it can be assumed they created a Q pattern (W/m3) which is the simulated ultrasound power propagation through the CT [10]. They showed success in a small number of patient models. When analyzing focal spot heating, they compared the hottest voxel of the MRTI and the simulation images. The hottest voxel is the single voxel around or at the focal spot where 10 temperature is the highest. The hottest voxel is important to determine how well the focusing worked as the higher the temperature, the more effective the phase correction focusing and the treatment. We will also utilize the hottest voxel method to look at temperature rise. In addition, we will model the acoustic properties of diploe in the skull in the study that follows and achieve results for all patient models. 11 METHODS The tcMRgFUS simulation was designed based on Insightec clinical treatments performed at the University of Utah. More treatments continue to take place, but for this study, data from 20 patients were used. The patient data were used under IRB approval and de-identified to protect patient privacy. Patients are referenced by ID number when different sets of specific patient-derived data are discussed. Clinical treatment data are used because they provided a direct way to match simulations to actual human treatments. These treatments were already performed so patient data were accessible and used to make the simulation. No patients were excluded from the study to test how capable the simulation was of matching the temperature rise in any actual skull model. Magnetic resonance images (MRI), CT, and treatment export files for each patient were all employed to create the simulations. The patient treatment export files contain necessary information for the simulation including transducer phases, phase steering, and power and energy output. This simulation takes place entirely in MATLAB [27] from inputting CT images used as the basis of the model to the thermal image outputs and data analysis. To begin, the CT scan of each patient’s head was first transformed into the transducer space. The transformation takes the CT image into the MRI orientation of the same patient. Transducer space and MR space are both the same as the patient lies in the transducer within the MRI scanner, therefore the coordinate systems are identical. Ensuring that the CT coordinates align with the MR guarantees that the focal spot in the CT image is the same, and that the acoustic simulation will function correctly. For the acoustic simulation to work, the CT is then segmented into seven different tissue types. Tissue types used are water, cortical bone, diploe, skin, gland, brain, and fat based on Hounsfield units. 12 Table 1 summarizes the acoustic values used in the acoustic simulation. The acoustic values used are for the speed of sound, attenuation, absorption, and density. The values are taken from the IT’IS Foundation acoustic value database [28]. These values were input into HAS for the seven different tissue types. The attenuation and absorption of diploe and bone were adjusted to best simulate expected results for the first patient modeled. Attenuation and absorption are the same for most values due to a scarcity of available information distinguishing between the two [28]. The CT segmentation, called the model, is at a resolution of 0.5 mm x 0.5 mm x 0.5 mm. A higher resolution allows a more accurate representation of heating as smaller voxel sizes provide more detail. Each segmented voxel in the segmented model is 0.5 mm3 and is represented by a different color. The model is assigned acoustic values for each of its segmented tissue types. An example of the segmented model can be seen in Figure 2. The threshold values used to create this segmentation can be seen in Table 2. The CT images were segmented into their 7 tissue types using the Hounsfield units listed within Table 2. Table 1: Acoustic properties used in HAS. Speed of sound (m/s) Attenuation (Np/cm/MHz) Absorption (Np/cm/MHz) Density (kg/m3) Water Gland Skin Fat Brain Diploe 1546.3 Cortical Bone 2653 1500 1560 1624 1440 0 0.133 0.22 0.044 0.068 1.3 1.3 0 0.133 0.22 0.044 0.012 1.3 1.3 1000 1064 1109 911 1046 1738 1300 1800 13 Figure 2: Patient 15 segmented CT. Axial, coronal, and sagittal views of the segmented CT at the focal point slice. At a 0.5 mm x 0.5 mm x 0.5 mm resolution and focal point is at the x, y, z location (226, 226, 300) in the center of each image. Table 2: Threshold values for the segmented models. Water Threshold <-200 (HU) Gland Skin Fat 50 >0 -200 to -100 to -101 Brain 299 Cortical Diploe Bone 1000+ 300 to 999 Acoustic Simulation Once the CT model is transformed and segmented, acoustic simulations take place. The acoustic simulation is performed in HAS. HAS uses the acoustic properties of the model to simulate and create a power deposition or Q pattern (W/m3). HAS was adapted to simulate a transcranial transducer by splitting up the transducer into seven zones as the transducer diameter is larger than the focal length. The HAS simulations are performed once for each of the three phase correction methods using each patient’s data. The phase corrections methods used are 1) no phase correction, 2) clinical phases, and 3) time reversal phase correction. HAS simulates no correction by setting the phase values to 0. Using clinical phases is performed by taking the phases used during the clinical treatments from the treatment export files and implementing them into HAS. If the CT model is in the 14 correct coordinate system with the appropriate acoustic values, this should match up best with the clinical treatments since they use identical phases. Time reversal [13] is a method of phase correction built into HAS which backward propagates from a focal point source back to the transducer to create the phases based on the model. The phases are based on the model and do not require data entry or input as with the clinical phase correction. Once a Q pattern is made for each phase correction method, these can be put into the thermal simulation to solve for temperature. Thermal Simulation The Q pattern is input into the thermal simulation to formulate tissue temperatures expected from the various models. The Pennes Bioheat equation [29] is used to create the temperatures by calculating temperatures based on the Q pattern and applying thermal conductivity, perfusion, diffusion, and specific heat of the skull. This is repeated individually for each method of phase correction. Once thermal images have been made, the images are transformed into the MRTI orientation and down-sampled into the MRTI resolution. This allows for the most accurate comparison between simulation and clinical images because the voxels are at the same size and orientation. MRTI images are taken in either the sagittal, coronal, or axial directions which also affects the amount of heating captured by the images. Image voxel size then becomes 1.094 mm x 2.188 mm x 3 mm after down-sampling. Simulation images can then be compared. 15 Statistical Analysis Simulations were compared primarily by the temperature increases of the hottest voxel. The hottest voxel was analyzed over time for each method of phase correction and for the MRTI, and also at the peak temperature. Each method of phase correction was compared in terms of temperature increase, and the MRTI was compared mainly with the clinical phase correction as both used the same phases. A regression analysis and t-test was performed between the MRTI and clinical phase correction to assess how closely the two compared for each patient, and to assess how patient-specific skull differences might affect this comparison. An analysis of variance (ANOVA) test was performed for the overall hottest voxel from all phase correction simulations normalized to the MRTI for the last two sonications of each patient. The energy of the transducer generally peaks during the last two sonications and the treatment concludes. The temperatures achieved then should be high enough to ablate the tissue at the focus. This tested for significant differences between the heating across phase correction methods. Once a significant difference was discovered, a Tukey honest significance difference (HSD) test was performed to observe which phase correction temperature rises were significantly different from one another. 16 RESULTS The results analyze how the implemented phase correction methods function in respect to the clinical MRTI images and to one another. They are analyzed over time and in terms of the temperature rise at the hottest voxel. A tighter focus generally provides a higher temperature rise as more energy is focused in one location. An example of the focus for the different phase correction methods can be seen in Figure 3. Time reversal appears to have the tightest focus while no phase correction is blurred and has a lower temperature rise. This hottest voxel is observed over time as well in Figures 4. The temperature rise in all three phase correction cases follows the same heating and cooling time frames as the MRTI. Figure 4 shows two sonications from 3 different patients that were simulated. Figures 4A and 4B depict an ideal simulation case where clinical phase correction and MRTI match closely during the entire heating and cooling. Figure 4C Figure 3: Hottest voxel thermal image comparison at the focus. 20 x 20 axial plane voxel view, with each voxel being .5 mm x .5 mm. Temperature rise is in terms of degrees Celsius. No correction shown heating the least with a blurred focal spot while time reversal had a tighter focus and produced the highest temperature rise. 17 A. B. C. D. E. F. Figure 4: Hottest voxel temperature rise over time. Temperature rise observed for no phase correction (No PC), clinical phase correction (Clinical PC), time reversal phase correction (TR PC), and MRTI. Figures 4A and 4B are from Patient 15’s sonication 17 and 18. Figures 4C and 4D are from Patient 17’s sonication 8 and 9. Figures 4E and 4F are of Patient 28’s sonications 4 and 5. MRTI is not a smooth line as it is made of 10-12 images over the acquisition duration. These 6 cases show how time reversal generally created the hottest temperature rise while no correction had the least. Clinical correction was between time reversal and no correction, and had a similar temperature rise to the MRTI. 18 contains a case where clinical correction peaked above time reversal, and Figure 4D has an MRTI that does not match the heating well. In order to validate the simulation further, the hottest voxel temperature rise for clinical phase correction and MRTI are analyzed further. These temperatures should match the best because they use the same exact phase aberration correction. If these phases match, the simulation would be successful in terms of replicating the clinical treatment with the current model. Figure 5 compared approximately half of the sonications for 6 different patients. The MRTI and clinical phase correction is compared for every sonications temperature rise MRTI time points. The MRTI has approximately 10-12 MRTI images per sonication, and the temperatures at the second through fourth images are compared to see how well the simulation matches the treatment. These time points always coincide with the temperature rise. The simulated phase correction methods contain approximately 10 times more images and are compared based on those three MRTI time points. The entire simulation is not compared as the heating is what is being primarily analyzed and the cooling of the simulation is not the focus at this time. The simulation would be unsuccessful if the first three images during the temperature rise failed to match. Table 3 summarizes the comparison of the simulated clinical temperature rise to the clinical MRTI. This compares the temperature rise of the MRTI to the simulated clinical correction using a t-test for all 20 patients and lists the R2 values. The result states if the test was significantly different. Four of the patients simulated were found to be statistically different in terms of matching the clinical MRTI to the simulated clinical phase correction. A t-test was performed for each patient to find any cases where the simulation may have not matched the clinical temperature rise. The R2 value for each patient was also found to 19 see how well the linear regression matched up between the simulated and clinical temperatures. In an ideal scenario, the R2 value would be 1 which would signify the simulated temperatures matching exactly to the clinical MRTI temperatures. A. B. C. D. E. F. Figure 5: Comparison of the hottest voxel temperatures of 6 different patients. The temperature comparison is between the MRTI and clinical phase correction for the first three time points of the last three hottest voxels. These 6 were selected randomly to show examples of high R2 values and low R2 values from the linear regressions for these time points. R2 line close to identity line in these 6 cases which show that MRTI temperature rise was similar to the clinical temperature rise. Figure 5A is from patient 15, 5B is from patient 17, 5C is from patient 21, 5D is from patient 23, 5E is from patient 27, and 5F is from patient 55. 20 Table 3: Comparison of simulated and clinical phase correction time points 2-4 for all patients. R2 Patient Total Sonications Slope P-value Test Result Number Sonications Used 15 18 10-18 0.9403 1.0560 0.8329 0.2121 Insignificant 17 11 6-11 0.9083 0.9318 0.6766 -0.4207 Insignificant 21 13 8-13 0.7813 0.8555 0.5589 -0.5903 Insignificant 23 10 5-10 0.9098 1.0836 0.9888 -0.0142 Insignificant 24 8 5-8 0.8775 0.8986 0.8540 -0.1852 Insignificant 27 12 6-12 0.8348 1.0346 0.6015 0.5264 Insignificant 28 10 5-10 0.8631 0.5368 0.0018 -3.3799 P<.05, significant 34 5 2-5 0.8698 1.4281 0.2076 1.2984 Insignificant 39 9 5-9 0.5957 0.5936 0.000025 5.0318 P<.05, significant 40 10 5-10 0.8133 1.1297 0.2907 1.0733 Insignificant 42 7 3-7 0.9008 2.0837 0.0194 2.4810 P<.05, significant 43 6 3-6 0.8467 1.2970 0.1417 1.5242 Insignificant 44 5 3-5 0.8404 1.8292 0.1148 1.6677 Insignificant 46 6 3-6 0.9247 1.0160 0.7429 0.3321 Insignificant 47 6 3-6 0.9815 2.6245 0.0138 2.6745 P<.05, significant 48 5 3-5 0.9487 1.3547 0.4514 0.7719 Insignificant 49 6 3-6 0.9134 1.1934 0.3398 0.9757 Insignificant 50 5 2-5 0.9452 2.1781 0.2391 1.2101 Insignificant 51 7 4-7 0.9508 0.9482 0.2631 1.1600 Insignificant 55 6 3-6 0.9101 0.9932 0.7441 -0.3305 Insignificant statistic 21 A slope of 1 would also signify if the temperatures match, which is reported in Table 3. A confidence interval of 95% was used for each t-test. 16 out of 20 t-tests determined that the temperature rise for 16 patients simulated was not statistically different from the MRTI. All methods of phase correction were lastly normalized to the MRTI and compared to one another. An Analysis of Variance (ANOVA) test was performed to compare the different methods of phase correction to one another to analyze if any of the means normalized to MRTI are statistically different. The ANOVA found at least one group was statistically different from another. A Tukey multiple comparison test was then performed between each phase correction method to determine if other pairings of results were statistically different from one another. Figure 6 shows the average hottest voxel for each method of phase correction normalized to MRTI. The average was calculated by finding the hottest voxel for the last two sonications of every patient and dividing each individual value by the hottest MRTI voxel for the same sonication. In the figure, no phase correction had the lowest average temperature rise, while time reversal had the highest. Clinical phase correction was in between time reversal and no correction. Table 6 summarizes the results of the Tukey HSD test performed which found all groups to be statistically different from one another. The temperature rise of time reversal compared to MRTI was on average 1.20 times higher, while no correction was only 0.68 as high as MRTI. Clinical correction was on average 0.97 that of the MRTI. 22 Figure 6: The average hottest voxel temperature rise normalized to MRTI. Hottest voxels taken from the last two sonications of all 20 patients phase correction simulations and averaged over MRTI. Time reversal temperatures are on average 120% higher than MRTI, while no phase correction was only 68% that of the MRTI. Clinical phase correction was similar in terms of temperature rise at 97% when normalized to MRTI. Table 6: Tukey HSD Test for temperature rise of all phase correction methods normalized to MRTI. Phase Correction Tukey HSD Q p-value Result Simulated Statistic Temperature Pair No PC vs Clinical 7.3793 0.0010053 ** p<0.01, significant No PC vs TR 13.3510 0.0010053 ** p<0.01, significant Clinical vs TR 5.9716 0.0010053 ** p<0.01, significant 23 DISCUSSION Current treatment methods of essential tremor are highly invasive and are risky for many elderly patients. Magnetic resonance-guided focused ultrasound (MRgFUS) provides an alternative treatment for Parkinson’s and essential tremor. MRgFUS is used to ablate neural tissue associated with tremor using acoustic energy. The acoustic energy is directed at the lesion from outside of the cranium. The sonication must then pass through both hard and soft tissues to safely eliminate target tissue. This becomes challenging due to energy defocusing caused by hard tissues. Defocusing can result in undesired heating of other tissues and ineffective treatments. Phase aberration correction provides the ability to improve ultrasound beam focusing in transcranial magnetic resonance-guided focused ultrasound (tcMRgFUS) treatments. By simulating treatments with phase correction, the outcome of a treatment can be visualized. The aim of this research was to simulate clinical tcMRgFUS treatments and apply different methods of phase correction to improve beam focusing. In addition, we aimed to match MRTI temperatures to clinical phase correction simulations to replicate clinical treatments. Transcranial MRgFUS phase aberration correction simulations can be validated by comparing thermal simulations of the different methods of phase correction to magnetic resonance thermal imaging (MRTI). Clinical phase correction and MRTI were found to be the most similar, indicating the simulation was successful at simulating the clinical conditions and replicating temperature rise. No phase correction was significantly worse in terms of the temperature rise compared to the other methods as expected, and time reversal correction had the highest temperature rise. Therefore, the simulation does well to match the clinical correction temperature rise to the real world MRTI. Also, clinical correction improves 24 beyond no phase correction’s temperature rise and time reversal provides higher simulated temperatures than both correction methods. Using phase aberration correction can improve the focusing of ultrasound therapies to provide a higher temperature rise. As expected, no phase correction proved to cause the least amount of temperature rise in comparison to the other methods of phase correction. No correction does not calculate phases and does not account for heterogeneities in the head physiology and anatomy that the ultrasound energy must pass through. It becomes aberrated by the geometry which can result in a significant loss of power at the focus. The average temperature rise simulated by no phase correction was only 68% as high as MRTI as seen in Figure 6. This was lower than the other two methods of phase correction. It simulated the worst and produced the lowest temperatures. Time reversal proved the opposite and produced higher temperatures. Time reversal phase correction calculated the highest temperature rise on average. Higher temperatures indicates that the phase correction is working and is providing a tight focal spot in the brain similar to that shown in Figure 3. The temperature rise of time reversal compared to MRTI was on average 120% higher. Time reversal phases were calculated based on the segmented model and should produce the highest temperatures. Clinical phase correction was calculated outside of HAS and was based on a different calculation method with a different original model. Therefore, time reversal should produce the highest temperatures in the simulation as confirmed in Figure 4 and Figure 6. Clinical correction was on average slightly lower than the MRTI with average temperatures 97% that of the MRTI. Time reversal and clinical correction simulations both had higher 25 amounts of error and the simulated correction was slightly lower than the measured MRTI. However, clinical phase correction and the MRTI were found to match well. No phase correction and time reversal were found to be different from clinical temperatures in the order expected. Interestingly, time reversal was not consistently higher than clinical correction. In Figure 4C, one instance of the temperature rise being reversed can be seen. Depending on the patient, time reversal had varying degrees of improvement from clinical correction. Some cases had time reversal being drastically higher than the clinical phases, while some others were close or reversed like in Figure 4C. Time reversal on average created the highest simulated temperatures followed by clinical correction and lastly, no phase correction. In current research, acoustic simulations have been performed for many different focused ultrasound applications. For tcMRgFUS, the most common methods are the Hybrid Angular Spectrum and k-Wave [11], [22] – [24], [30]. Both methods produce a pressure pattern from ultrasound beams propagating through a model and can be performed in MATLAB. Both methods can be used in transcranial applications. One example of kWave being used is in the McDannold paper [21]. Similar methods were used in the McDannold paper but with a different end goal [21]. This paper looked to form thermal simulations of tcMRgFUS but was focused on heating in the skull bone rather than in the focal spot in the brain. They also used k-Wave for the acoustic simulation instead of using HAS [11], [30]. The thermal simulation performed in this research could potentially be used to look at the damage in the skull. It would be interesting to see if there is any differences between these two acoustic simulations and the results they produce. 26 The Leung paper used HAS in their acoustic simulations and was the most similar to this research [26]. While the methods used were the same down to the software used, the Leung paper only showed a small number of patients. They provided results on a handful of patient models where this research has used 20 at this time and plans to use more. Both good and bad results were used while they only showed their most successful examples. By showing all patients simulated, an estimate as to the strength of the simulation accuracy could be observed when using the same methods as Leung, The model used has several limitations that could be further improved. The model is not an exact model of the head as it is a simple segmentation and excludes media like gray brain matter, white brain matter, and cerebral spinal fluid. The model also groups similar tissues together and assigns them singular acoustic values when they are far more complex and would have a range of acoustic values in each tissue type. In addition to the segmentation not being exact, the attenuation values for cortical bone and diploe were adjusted to accurately fit the first patient simulated. The attenuation values need to be further validated so that better values can be incorporated. More patients need to be run as well to obtain a larger sample size. At the time of this research, 20 patients were completed. There are approximately 60 total patient data sets available to be simulated. More patients could change the results and may require further adjustment of the parameters used in the simulation. The simulation also needs to be sped up as currently it takes 6 hours for a patient with 11 sonications to be run. For patient 15 who had 18 total sonications, the simulation took closer to 8 hours. This makes it difficult to quickly run the simulation to see treatment outcomes or to generate time reversal phase correction phases. 27 Improving the speed at which the simulation is run would make the software more suitable in treatment planning. Being able to run a treatment based off a computed tomography (CT) scan would make it possible to see how successful tcMRgFUS may be for a patient. Problems could be solved before even running a treatment. The simulation could also improve treatment phases and improve the overall accuracy of the treatment. A couple limitations in the simulation structure need to be improved upon before running for patients. Several new steps can be taken to improve the simulation and the phase correction methods. The parameters and model involved can be improved further through more accurate acoustic values and more detailed models. The acoustic values of bone need to be validated with hydrophone scans of ex vivo skull fragments [31], [32]. Hydrophones measure the acoustic field produced by a focused ultrasound transducer. Hydrophone scans can be used to measure the acoustic values of bone and those values would replace the current values in the simulation. Inaccuracies in these values can skew the predicted temperatures [33]. The current attenuation and absorption values used for bone and diploe were estimated based on the first patient and need to be further validated in the future. More clinical treatment patients will also be run through the simulation and the parameters can be further adjusted once a larger data set is accumulated. A different way to model the patients will be implemented, making the models more detailed. Rather than segmenting the skull and then assigning acoustic properties, other studies have assigned speed of sound and attenuation as continuous functions of the CT HU [26]. No segmentations would need to be performed and the CT of the patient would be able to be simulated based on the exact structure. This could make tcMRgFUS further specific to the patient. 28 The research conducted can potentially improve focusing and allow for more patient-specific treatments. This may improve treatment specificity and help to spare surrounding brain tissue from unnecessary damage. Being able to simulate potential treatments can assist with treatment planning and efficiency to view a treatment before implementation. Furthermore, using these simulations to find the best phase correction method helps ensure each treatment will be maximally effective. Improved phase correction methods where more energy reaches the focus have the potential to allow for patients with low skull density ratios to be included in tcMRgFUS treatments. 29 REFERENCES [1] "Tremor Fact Sheet", NINDS, Apr. 2021, NIH Publication No. 21-NS-4734 [2] J. Jankovic, “Distinguishing Essential Tremor From Parkinson’s Disease”. Practical Neurology, pp. 36-38, 2012. [3] “Essential tremor”, Mayo Clinic, MFMER, Jan. 2019. [4] N. Y. Jung and J. W. Chang, “Magnetic resonance-guided focused ultrasound in neurosurgery: Taking lessons from the past to inform the future,” J. Korean Med. Sci., vol. 33, Oct. 2018. [5] V. Krishna, F. Sammartino, and A. Rezai, “A review of the current therapies, challenges, and future directions of transcranial focused ultrasound technology advances in diagnosis and treatment,” JAMA Neurology, vol. 74, pp. 246-254, Feb. 2018. [6] P. Ghanouni et al., “Transcranial MRI-guided focused ultrasound: A review of the technologic and neurologic applications,” Am. J. Roentgenol., vol. 205, pp. 150159, Jul. 2015. [7] A. Boutet et al., “The relevance of skull density ratio in selecting candidates for transcranial MR-guided focused ultrasound,” J. Neurosurg., vol. 132, pp. 17851791, May 2019. [8] F. J. Fry and J. E. Barger, “Acoustical properties of the human skull,” J. Acoust. Soc. Am., vol. 63, pp. 1576-1590, May 1978. [9] S. Almquist, D. L. Parker, and D. A. Christensen, “Rapid full-wave phase aberration correction method for transcranial high-intensity focused ultrasound therapies,” J. Ther. Ultrasound, vol. 4, Dec. 2016. 30 [10] S. Almquist, J. de Bever, R. Merrill, D. Parker, and D. Christensen, “A full-wave phase aberration correction method for transcranial high-intensity focused ultrasound brain therapies,” in Annu Int Conf IEEE Eng Med Biol Soc, 2014, pp. 310-313. [11] B. E. Treeby, J. Jaros, A. P. Rendell, and B. T. Cox, “Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method,” J. Acoust. Soc. Amer., vol. 131, pp. 4324-4336, Jun. 2012. [12] M. Karaman, A. Atalar, H. Koymen and M. O'Donnell, "A phase aberration correction method for ultrasound imaging," in IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 40, pp. 275-282, July 1993. [13] M. Fink, “Time Reversal of Ultrasonic Fields—Part I: Basic Principles,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 39, pp. 555-566, Sep. 1992. [14] T. J. Hudson et al., “Simulating thermal effects of MR-guided focused ultrasound in cortical bone and its surrounding tissue,” Med. Phys., vol. 45, pp. 506-519, Feb. 2018. [15] N. Sharma and L. M. Aggarwal, “Automated Medical Image Segmentation Techniques,” Med. Phys, vol. 35, pp. 3-14, Jan. 2010. [16] F. Gaillard and K. Greenway, “Hounsfield unit.” Radiopaedia, https://radiopaedia.org/articles/hounsfield-unit?lang=us [17] T. D. Webb et al., “Measurements of the Relationship between CT Hounsfield Units and Acoustic Velocity and How It Changes with Photon Energy and 31 Reconstruction Method,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 65, pp. 1111-1124, Jul. 2018. [18] J. Aarnio, G. T. Clement, and K. Hynynen, “A new ultrasound method for determining the acoustic phase shifts caused by the skull bone,” Ultrasound Med. Biol., vol. 31, pp. 771-780, Jun. 2005. [19] M. Tanter, J.-L. Thomas, and M. Fink, “Focusing and steering through absorbing and aberrating layers: Application to ultrasonic propagation through the skull,” J. Acoust. Soc. Am., vol. 103, pp. 2403-2410, May 1998. [20] Y. G. Kim, E. J. Kweon, W. S. Chang, H. H. Jung, and J. W. Chang, “Magnetic Resonance-Guided High Intensity Focused Ultrasound for Treating Movement Disorders,” Prog Neurol Surg, vol. 33, pp. 120-134, Jan. 2018. [21] N. McDannold, P. J. White, and R. Cosgrove, “Predicting bone marrow damage in the skull after clinical transcranial MRI-guided focused ultrasound with acoustic and thermal simulations,” IEEE Trans. Med. Imaging, vol. 39, pp. 3231-3239, Apr. 2020. [22] U. Vyas and D. Christensen, “Ultrasound beam propagation using the hybrid angular spectrum method,” in Annu Int Conf IEEE Eng Med Biol Soc, 2008. [23] U. Vyas and D. Christensen, “Extension of the angular spectrum method to calculate pressure from a spherically curved acoustic source,” J. Acoust. Soc. Am., vol. 130, pp. 2687-2693, Nov. 2011. [24] U. Vyas and D. Christensen, “Ultrasound beam simulations in inhomogeneous tissue geometries using the hybrid angular spectrum method,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, vol. 59, pp. 1093-1100, Jun. 2012. 32 [25] U. Vyas, P. Ghanouni, C. H. Halpern, J. Elias, and K. B. Pauly, “Predicting variation in subject thermal response during transcranial magnetic resonance guided focused ultrasound surgery: Comparison in seventeen subject datasets,” Med. Phys., vol. 43, pp. 5170, Sep. 2016. [26] S. A. Leung, T. D. Webb, R. R. Bitton, P. Ghanouni, and K. Butts Pauly, “A rapid beam simulation framework for transcranial focused ultrasound,” Sci. Rep., vol. 9, pp. 7965, 2019. [27] MATLAB. (2020). version 9.9.0 (R2020b). Natick, Massachusetts: The MathWorks Inc. [28] Hasgall PA, Di Gennaro F, Baumgartner C, Neufeld E, Lloyd B, Gosselin MC, Payne D, Klingenböck A, Kuster N, “IT’IS Database for thermal and electromagnetic parameters of biological tissues,” Version 4.0, May 15, 2018, DOI: 10.13099/VIP21000-04-0. itis.swiss/database [29] H. H. Pennes, “Analysis of tissue and arterial blood temperatures in the resting human forearm,” J. Appl. Physiol., vol. 1, pp. 93-122, Aug. 1948. [30] J. Jaros, F. Vaverka, and B. E. Treeby, “Spectral domain decompo-sition using local Fourier basis: Application to ultrasound simulation on a cluster of GPUs,” Supercomput. Frontiers Innov., vol. 3, pp. 40–55, Nov. 2016. [31] G. Pinton, J. F. Aubry, E. Bossy, M. Muller, M. Pernot, and M. Tanter, “Attenuation, scattering, and absorption of ultrasound in the skull bone,” Med. Phys., vol. 39, pp. 299-307, Jan. 2012. 33 [32] S. Pichardo, V. W. Sin, and K. Hynynen, “Multi-frequency characterization of the speed of sound and attenuation coefficient for longitudinal transmission of freshly excised human skulls,” Phys. Med. Biol., vol. 56, pp. 219-250, Jan. 2011. [33] T. E. Vaughan and K. Hynynen, “Effects of parameter errors in the simulation of transcranial focused ultrasound,” Phys. Med. Biol., vol. 47, pp. 37-45, Jan. 2002. 34 Name of Candidate: Emma N. Slominski Birth date: December 18, 1998 Birth place: Rockville, Maryland Address: 470 S 1300 E Salt Lake City, UT 84102