Deep Brain Ultrasound Ablation Thermal Dose Modeling with in Vivo Experimental Validation
Zhanyue Zhao Benjamin Szewczyk Matthew Tarasek Charles Bales Yang Wang Ming Liu Yiwei Jiang Chitresh Bhushan Eric Fiveland Zahabiya Campwala Rachel Trowbridge Phillip M. Johansen Zachary Olmsted Goutam Ghoshal Tamas Heffter Katie Gandomi Farid Tavakkolmoghaddam Christopher Nycz Erin Jeannotte Shweta Mane Julia Nalwalk E. Clif Burdette Jiang Qian Desmond Yeo Julie Pilitsis Gregory S. Fischer






Journal of  Class Files, Vol. 14, No. 8, August 2021
Shell et al.: A Sample Article Using IEEEtran.cls for IEEE Journals

Deep Brain Ultrasound Ablation Thermal Dose Modeling with in Vivo Experimental Validation
    Zhanyue Zhao, Benjamin Szewczyk, Matthew Tarasek, Charles Bales, Yang Wang, Ming Liu, Yiwei Jiang, Chitresh Bhushan, Eric Fiveland, Zahabiya Campwala, Rachel Trowbridge, Phillip M. Johansen, Zachary Olmsted, Goutam Ghoshal, Tamas Heffter, Katie Gandomi, Farid Tavakkolmoghaddam, Christopher Nycz, Erin Jeannotte, Shweta Mane, Julia Nalwalk, E. Clif Burdette, Jiang Qian, Desmond Yeo, Julie Pilitsis, and Gregory S. Fischer

Z. Zhao, B. Szewczyk, C. Bales, Y. Wang, M. Liu, Y. Jiang, K. Gandome, F. Tavakkolmoghaddam, C. Nycz, and G. S. Fischer are with Worcester Polytechnic Institute, Worcester, MA e-mail: zzhao4@wpi.edu, gfischer@wpi.edu.

B. Szewczyk, J. Qian, and J. Pilitsis are with the Department of Neurosurgery, Albany Medical Center, Albany, NY

M. Tarasek, C. Bhushan, E. Fiveland, and D. Yeo are with GE Global Research Center, Niskayuna, NY

Z. Campwala, R. Trowbridge, Z. Olmsted, S. Mane, J. Nalwalk, J. Qian, and J. Pilitsis are with the Department of Neuroscience and Experimental Therapeutics, Albany Medical Center, Albany, NY

P. M. Johansen and J. Pilitsis are with Charles E. Schmidt College of Medicine, Florida Atlantic University, Boca Raton, FL

E. Jeannotte is with Animal Resources Facility, Albany Medical Center, Albany, NY

G. Ghoshal, T. Heffter, and E. C. Burdette are with Acoustic MedSystems, Inc., Savoy, IL

This research is supported by National Institute of Health (NIH) under the National Cancer Institute (NCI) under Grant R01CA166379 and R01EB030539.

    
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



§ ABSTRACT 

Intracorporeal needle-based therapeutic ultrasound (NBTU) is a minimally invasive option for intervening in malignant brain tumors, commonly used in thermal ablation procedures. This technique is suitable for both primary and metastatic cancers, utilizing a high-frequency alternating electric field (up to 10 MHz) to excite a piezoelectric transducer. The resulting rapid deformation of the transducer produces an acoustic wave that propagates through tissue, leading to localized high-temperature heating at the target tumor site and inducing rapid cell death. To optimize the design of NBTU transducers for thermal dose delivery during treatment, numerical modeling of the acoustic pressure field generated by the deforming piezoelectric transducer is frequently employed. The bioheat transfer process generated by the input pressure field is used to track the thermal propagation of the applicator over time. Magnetic resonance thermal imaging (MRTI) can be used to experimentally validate these models. Validation results using MRTI demonstrated the feasibility of this model, showing a consistent thermal propagation pattern. However, a thermal damage isodose map is more advantageous for evaluating therapeutic efficacy. To achieve a more accurate simulation based on the actual brain tissue environment, a new finite element method (FEM) simulation with enhanced damage evaluation capabilities was conducted. The results showed that the highest temperature and ablated volume differed between experimental and simulation results by 2.1884^∘C (3.71%) and 0.0631 cm^3 (5.74%), respectively. The lowest Pearson correlation coefficient (PCC) for peak temperature was 0.7117, and the lowest Dice coefficient for the ablated area was 0.7021, indicating a good agreement in accuracy between simulation and experiment.




Needle-Based Therapeutic Ultrasound (NBTU), Magnetic Resonance-Guided Robotically Conformal Ablation, Finite element modeling, in vivo Swine Ablation




§ INTRODUCTION

Thermal ablation techniques are being developed as minimally invasive or non-invasive alternatives to primary and metastatic brain tumor surgery <cit.>. Different types of technologies were developed, including four conventional methods: (1) Radio frequency (RF) ablation, which uses electricity from an alternating current at a favorite frequency of around 500kHz to produce ionic agitation in tissue, which results in friction and heating of tumor <cit.>. (2) Microwave ablation (MWA), which uses electromagnetic waves up to 2450MHz to generate polar molecules agitation with friction and heat in tissue and causes coagulative necrosis of cells <cit.>. (3) Laser ablation, which utilizes intense and highly-collimated beams of monochromatic light that are designed for biological tissues for thermal destruction of tumors <cit.>. (4) Ultrasound ablation, which uses the transmission of high-frequency sound waves through biological tissue transfers mechanical energy that leads to localized heating and coagulative necrosis of surrounding cells <cit.>. The ultrasound ablation contains two types of techniques, namely extracorporeal methods using high intensity focused ultrasound (HIFU) for non-invasive ablation procedure <cit.>, and intracorporeal method using an inserted transducer with direct contact of the target area for minimally invasive ablation procedure <cit.>. There is a promising new ablation technology developed in recent years called pulsed-field ablation (PFA), which uses a train of microsecond duration high amplitude electrical pulsed that ablates myocardium by electroporation of the sarcolemmal membrane without measurable tissue heating, and this method is aiming for the treatment of atrial fibrillation<cit.>.

Numerical modeling of thermal ablation solves the physical equations governing energy deposition and heat transfer in tissue to determine the transient temperature profile and assess tissue damage upon heating <cit.>. However, due to patient-specific tissue compositions that might affect local features like absorption, attenuation, and perfusion, thermal ablation can be challenging to anticipate with absolute accuracy. Finite element modeling (FEM) can be applied in thermal ablation modeling. It can calculate the piezoelectric transducer's acoustic pressure and the transferring bioheat deposition along the desired acoustic medium by subdividing a large system into smaller ones. It can also simpler elements by meshing the construction. The finite element method formulation of a boundary value problem finally results in a system of algebraic equations. 

In this work, an extended FEM simulation towards NBTU application in a real brain tissue environment with tissue damage evaluation is one of the primary studies for MR-guided robot-assisted conformal ablation procedures. We developed a Comsol FEM simulation based on a real brain tissue medium with a CEM43 isodose map for evaluating tissue damage. Moreover, we validated the model with in vivo swine ablation experimental results and proved the feasibility and accuracy of this new model.



§ FEM ANALYSIS OF NBTU IN BRAIN TISSUE

Temperature measurements can provide absolute or relative temperature changes in comparison to unheated tissue or reference data<cit.>. However, thermal damage is generally evaluated using the Sapareto-Dewey equation, and the ablation zone is defined by the cumulative equivalent minute (CEM) standard<cit.>. The CEM43 is defined as the equivalent time the target tissue temperature has to be at 43^∘C to induce thermal damage<cit.>. We selected a CEM43 of 70 based on our collaborators' experience with ablation in liver tissue<cit.>,  also a comparison to similar paper<cit.> and studies demonstrating a CEM43 of 50 to 240 range to induce cell death with soft tissue<cit.>. Using Sapareto and Dewey's formulation, the thermal damage dosage was computed using the temperature from the MRTI given by Equation <ref>


    CEM43=∑_t=0^t=final R^43-T_tΔ t, {    R=0.25 for T<43^∘ C
        R=0.50 for T≥43^∘ C
            .


where T_t is the average temperature during time Δ t. The unit of thermal dose is equivalent to minutes at 43^∘C. The application of CEM43 can estimate the damage to brain tissue.

Based on the thermal damage of soft tissue, a FEM simulation can be created and used as a numerical prediction tool for accurately tracking thermal changes using NBTU in preparation for human trials.



 §.§ FEM Analysis of NBTU in Brain Tissue


The FEM study for simulating thermal propagation and CEM43 threshold calculation was informed by NBTU parameters during surgery, while the parameters for the application were based on the phantom study <cit.>. We enhanced our existing numerical models with an acoustic medium based on brain-tissue-specific parameters for comparison to our in vivo animal data. 



  §.§.§ Simulation Setup


The FEM model was implemented in Comsol to model the probe’s mechanical deformation, resultant applicator stationary acoustic pressure field, and the thermal damage isodose threshold of 70 CEM43. A cross-sectional model of the probe’s geometry was created with a dimensional parameter of 1.5mm outer diameter by 1.1mm inner diameter ring, and our notches of 0.1mm depth by 0.05mm width were also added such that the probe was segmented into 90^∘ and 180^∘ sectors. A detailed drawing of the probe dimension is shown in Fig. <ref>. Lead zirconate titanate (PZT-4) from Comsol material library was selected as the probe cylinder material with no further modifications and the simulation properties for the transducer were used as shown in Table <ref>. The probe was surrounded by a 100 mm × 100 mm (L×W) 2D acoustic medium. 





The solid mechanics physics interface defined the transducer as a piezoelectric material and defined a fixed constraint within the inner surface of the cylinder. The piezoelectric poling direction was defined by a base vector coordinate system defined by the transform shown in Table <ref> and was aligned radially outward from the center of the transducer geometry. The electrostatics physics interface was also used to define the electric potential, or electrical terminal along the 90^∘ outer surface of the transducer and to define a ground constraint along the inner surface. Different types of probes (ablation range) were used in the experiment, and the specific parameters of probes that were used in the simulation are shown in Table <ref>. The electric potential can be calculated by Equation <ref>:


    V = √(P(1+ε)Z/η)


Where P is the desired acoustic power we need for ablation, ε is the power loss generated from the RF inline filter, and the power loss for these probes has a range of approximately 10-20% based on the probe type and the method of fabrication, based our collaborator's suggestion we used median value of 15% in our simulation. Z is the transducer impedance value, and we used 50ohm as a general piezoelectric material value. η is the efficiency of the probe, detailed value can be found in Table <ref>.





The pressure acoustics physics interface was used to define the wave propagation through the medium as a linear elastic fluid model with an attenuation coefficient of 31.96 Np/m. The bioheat transfer physics interface was used to define a heat source applied to the medium with a user-defined activation function shown in equation 1, where Q is the absorbed ultrasound energy over the acoustic field, acpr. Q_pw is the dissipated power intensity, step is a step function from 0 to 1 with smoothing of 0.005, t is the current time step of the simulation, and t_ProbeOn is a constant that represents the total time the probe was activated.


The brain tissue referenced acoustic medium material properties are described in Table <ref>. Far-field and thermally insulated boundary conditions were applied to the acoustic medium, and no reflection of ultrasonic waves from edges was assumed. A frequency domain study was conducted to obtain the acoustic pressure field produced by the applicator, applying a tetrahydral mesh with a maximum element size of λ/6, where λ represented the wavelength of the produced ultrasonic wave within the medium.





  §.§.§ Acoustic Pressure Mapping


The generated ultrasound waves produced by the transducer can be simulated in Comsol using a frequency domain study at the selected resonant mode, and the pattern mapping can be found in Fig. <ref>. Considering the properties of brain tissue performed small change in speed of sound propagation in both mediums, the patterns are similar to our previous simulation work reported in <cit.>.






  §.§.§ CEM43 Mapping



The bioheat transfer time-dependent study was conducted to validate the heat deposition by the simulated acoustic pressure field, and we modified the default equation in this node by using the generalized Penne’s bioheat transfer equation shown in Equation <ref>:


    ρ C∂ T/∂ t=∇· k ∇ T-c_bρ_bω_b(T-T_b)+Q_m+Q_p


Where ρ (kg/m^3) is the tissue density, C (J/kg/^∘ C) is the tissue specific heat capacity, k (W/m/^∘ C) is the tissue thermal conductivity, T (^∘ C) is the current temperature, T_b (^∘ C) is the temperature of blood, and we assumed the initial temperature of both tissue and blood as 37^∘C. ρ_b (kg/m^3) is the blood density, ω_b (kg/m^3/s) is the blood perfusion, c_b (J/kg/^∘ C) is the blood specific heat capacity, which we used the same value as the brain tissue heat capacity (C). Q_p (W/m^3) is the heat deposition due to sonication given by Equation <ref>


    Q_p=α_attenp^2/ρ c


Where α_atten (Np/m) is the medium attenuation, p (Pa) is the acoustic pressure field, ρ (kg/m^3) is the medium density, c (m/s) is the speed of sound within the medium. Q_m (W/m^3) is metabolic heat, which is the heat generated by the tissue. 

A general form PDE study was then conducted to calculate the CEM43 thermal dose map and 70 CEM43 isodose lines, which are shown in Fig. <ref>. Predicted thermal dose maps for NBTU probes showed inhomogeneous CEM 43 thermal dose patterns because of the random blood perfusion and metabolic heat influence generated by Comsol. Using the same volume acquisition method for MRT images, we used surface integration to calculate the 70 CEM43 isodose lines surrounding the area of the transducer in the middle section. Areas of additional slides with 5 mm slice thickness were acquired, multiplied by 5 mm, and summed together to calculate the numerical predicted ablated volume. 





  §.§.§ Model Sensitivity Analysis


We also performed a Taguchi-based sensitivity analysis with the above model to determine the contribution and ranking of parameters. We considered the maximum obtained temperature as the response variable, and the six 3-level parameters, namely Thermal Conductivity (k),  Heat Capacity (C), Metabolic Heat (Q_m), Blood Perfusion (ω_b), Attenuation (α_atten), and Voltage (V) applied onto the transducer, as factors of interest. The 3-level of each parameter based on human soft tissue proprieties, namely low, medium, and high, were shown in Table <ref>. To obtain the minimum number of simulations needed for the parameters with corresponding levels we used Taguchi analysis standard L-27 orthogonal array. After obtaining all combinations of parameters and their corresponding maximum temperature, we then performed the analysis of variance (ANOVA) to quantify the effect of each factor on the response variable <cit.>. Calculations of Taguchi design and ANOVA were all carried out in Minitab 21 software (Minitab LLC., PA). Results of ANOVA are shown in Table <ref>. As a result, V majorly affects the maximum obtained temperature with more than 50% contribution in the model. The second and third major contributing factors are α_atten and k, respectively. The aforementioned three parameters contribute over 83% in total to the response variable. C, Q_m, and ω_b are not significant in our model. 












§ MRTI IMAGING VALIDATION


This section summarized the experiments from the year 2020 to 2022 and selected swine data was used for validation. All the subject was approved by the Albany Medical College (AMC) Institutional Animal Care and Use Committee (IACUC). MR delivery of NBTU was conducted on male and female Sus scrofa domestics swine (8 to 20 weeks old, 10 weeks on average) weighing between 18 and 25kg. One acute swine (sacrificed immediately post-procedure), one subacute swine (sacrificed 1 to 4 days post-procedure), and 5 survival swines were selected to create a total of 7 lesions. MRT-images acquired from a 3T GE Signa Architect scanner (GE Global Research, Niskayuna, NY, USA) and the ablated volume calculated from MRT-images were compared with simulation results <cit.>. To operate the interstitial thermal ablation procedure, a seven-DOF NeuroAblation robot was created as an MR-compatible device. The end effector of the NeuroAblation robot is modular and may be changed to accommodate the preferred applicator for the intended intervention. In this work, a needle-based therapeutic ultrasound (NBTU) probe developed by Acoustic MedSystems Inc (Illinois, United States) powered by the TheraVision control system was deployed in the subject. A detailed description can be found in our previous work <cit.>.



 §.§ In-Vivo Swine Ablation Experiment



In this study, pigs were positioned in a prone position using a custom head holder to align with a NeuroAblation robot, which was fixed to MRI bed rails and registered with the MRI scanner's coordinate system. The procedure involved dividing the MRI room into a sterilized zone for the surgeon and nurse and a standard zone for the robot operator. The entry and target points for the procedure were selected using 3D MRI images, and the NeuroAblation robot validated the coordinates.

Once the points were confirmed within the robot's range, the pig was removed from the scanner, and the robot was operated to reach the target position, except for needle insertion, which was performed by the surgeon. Afterward, the pig was repositioned in the scanner for the ablation procedure. Thermal dose mapping was done using intraoperative MRI to track temperature changes and assess the ablation's effectiveness.

The procedure included low-dose test ablations to confirm the ablation volume's shape and direction before full treatment. Post-ablation imaging was conducted, and pigs were euthanized at various time points for further analysis. The study measured temperature distribution, thermal damage, and ablated volumes to assess the effectiveness of the procedure. A detailed description can be found in <cit.>.



 §.§ MRT-Images Temperature Validation


7 swine experimental data were selected for simulation validation. Firstly, the 90^∘ probe with 6w acoustic power under 690s duration time of ablation (swine 7 experimental and simulation results) was selected as one of the examples of analysis in this section, which consisted of the short period of low power ablation for pre-testing, and a high power ablation range, and the MRTI peak temperature change over time and compare with the simulation maximum temperature results are shown in Fig. <ref>. Note that the ablation process concludes both heating up and cooling down duration because during the cooling down period, the surrounding area was still being heating up, and the heating continued contributing to the CEM43 thermal isodose map. Results indicate that the center ablation surface performed between slice 2 and 3 considering only the two slices changed temperature significantly while the other slices changed very small, and the maximum temperature increased up to 57^∘C. The right figure shows the selected slices (slice 2 and 3) versus the simulation result in temperature change over time. Based on the simulation the predicted maximum temperature reached up to 58.5823^∘C, which yields maximum peak temperature deference with 1.5823^∘C (2.78%). The Pearson correlation coefficient (PCC) metrics were evaluated to compare the experimental and simulation results. The PCC between experimental data from slice 2 and the simulation peak temperature curve reached 0.9527.



The rest of the lesions' results are also compared with the simulation data, and all the images can be found in Fig. <ref> and Fig. <ref>. The results of 180^∘ is shown in Fig. <ref>, where the (a)-(b) shows the Swine 4 results with 180^∘ probe under 3W acoustic power and 180s duration time condition, and the peak temperature reached 62^∘C in experiment versus 63.2925^∘C in simulation with Pearson correlation coefficient (PCC) of 0.9566. (c)-(d) shows the Swine 6 results with 180^∘ probe under 4W acoustic power and 180s duration time condition, and the highest temperature reached 64^∘C in experiment versus 64.9839^∘C in simulation with PCC of 0.9275. Finally, the results of 360^∘ are shown in Fig. <ref>, where the (a)-(b) shows the Swine 1 experimental and simulation results with 360^∘ probe under 3W acoustic power and 90s duration time condition, and the highest temperature reached 41^∘C in experiment versus 42.3701^∘C in simulation with PCC of 0.9652. Note that the base temperature was 38^∘C. (c)-(d) shows Swine 2 results with 360^∘ probe under 3W acoustic power and 100s duration time condition, and the highest temperature reached 43^∘C in experiment versus 44.3343^∘C in simulation  with PCC of 0.9321. Note that there was no low-temperature testing process in this experiment. (e)-(f) Swine 3 results with 360^∘ probe under 3W acoustic power and 120s duration time condition, and the highest temperature reached 61^∘C in experiment versus 60.5230^∘C in simulation with PCC of 0.8145. Note that the base temperature was also 38^∘C in this experiment. (g)-(h) Swine 5 results with 360^∘ probe under 4W acoustic power and 120s duration time condition, and the highest temperature reached 59^∘C in experiment versus 61.1884^∘C in simulation with PCC of 0.7117. Detailed results can be found in Table <ref>. 








 §.§ Ablated CEM43 Map Results


The previous section discussed the peak temperature change between MRTI data and simulation results, however, the CEM43 isodose is the desired indicator for evaluating the tissue damage in the field <cit.>. By using the volumetric MRTI measured the temperature changes at different distances from the focal point, and the necrotic zone boundaries were measured at CEM4370 and multiplied the thickness of 5mm, we can calculate the ablated volume to compare with our simulation results.

The CEM43 of 70 isodose maps from 7 swine experiments are compared with the simulation results. An area was drawn around the margins of the ablated area with the lasso feature and then the area of the circled region in cm^2. Each ablated area was multiplied by the 5mm slice thickness if the ablation went through the posterior side of the slice. In slice 3 the largest distance from the probe center to the edge of the ablated area is 9mm and the ablated volume is 0.2344 cm^3, while in slice 3 the largest distance from the probe center to the edge of the ablated area is 7.4mm and the ablated volume is 0.2234 cm^3. Fig. <ref> shows the simulation versus the experimental results of the CEM43 of 70 isodose maps. The largest distance from the center of the probe to the edge of the ablated area is 9mm in the experiment versus 8.9778 in the simulation, which yields within 0.25% error. Another metric, the Dice coefficient was used for evaluating the degree of overlap between simulation and experimental CEM43 isodose pattern. For swine 7, the Dice coefficient between simulation (center slice of transducer) and experimental results (MTRI slice 2 data) is 0.8496. 



Fig. <ref> and Fig. <ref> show the simulation versus the experimental results of CEM43 70 isodose maps using 360^∘ and 180^∘ probe results, and the Dice coefficient of selected swines experiment versus the simulation can be found in Table <ref>. The lowest Dice coefficient of experiment and simulation results by using 360^∘ and 180^∘ are 0.7021 and 0.8145 respectively. Note that these results contain all the in vivo animal experiments the author had attended, which is beyond the 7 lesions mentioned in this study.
 




Detailed results of ablated volume data from the 7 swine experiments can be found in Table <ref>. By using the same method mentioned in MRTI, the ablated area calculated from simulated CEM43 of 70 isodose maps was multiplied by the 5mm slice thickness and summed together to get the ablated volume, which is also shown in Table <ref>. 





§ CONCLUSION


In this work, we developed a 2D version simulation with real brain tissue parameters focused on thermal damage modeling, and validated with multiple in vivo animal surgery. 

Taguchi-based sensitivity analysis was performed and was concluded that the voltage applied (V) and attenuation (α_atten) are the most important parameters to determine the obtained maximum temperature in our bioheat transfer model. Consistent with other researchers' findings, we observed significant performance of “voltage" contribution to ablation procedure <cit.>. We also suggested prioritizing attenuation (α_atten) and thermal conductivity (k) and carefully considering these parameters in both simulation and experiments in bioheat transfer. 

The 2D thermal damage FEM simulation in brain tissue, including the validation with multiple in vivo animal surgery. This model focused on real tissue thermal damage (CEM) instead of thermal propagation, which is what we desired in NBTU thermal ablation therapy. The highest temperature and ablated volume largest difference error between experimental and simulation results yield 2.1884^∘C (3.71%) and 0.0631cm^3 (5.74%) respectively, with the lowest Pearson correlation coefficient (PCC) of peak temperature is 0.7117, and the lowest Dice coefficient of ablated area is 0.7021. These results show good agreement of accuracy between simulation and experiment. Note that the ones with higher temperature error (swine 5, 3.71%) are not the same as the ones with higher volume error (swine 6, 5.74%), this may be because the temperature error is more dependent on thermal conductivity and absorption effect, while the volume error is more dependent on thermal dissipation and perfusion rate.




§ DISCUSSION


In this study, the changes of tissue properties are not considered, however as reported in <cit.> the perfusion and attenuation will change dynamically during the ablation procedure, especially the ablated section will be different compared to un-ablated status both in healthy tissue and tumor cells. Both low PCC and Dice coefficient were performed when using 360^∘, especially the heating up duration and some unexpected edge of ablation. Considering the sensitivity analysis we performed, the attenuation is ranking 2^nd among parameters in our simulation, and it can change up to ± 50% during the ablation procedure in soft tissue <cit.>, which will affect the ablation performance significantly. Also, the surrounding vasculature system may further induce these unexpected misalignments of ablation boundaries. Moreover, the heat transfer method between the above two different status tissues should also be considered to make the modeling more accurate. Future work may focus on dynamic parameters integrated simulation, especially considering the value shift from attenuation, and add more physics nodes with heat transfer among blood vascular system. The current sensitivity study only focuses on multiple parameters. In the future, we may consider involving some interactions of parameters in the analysis to get a more comprehensive ranking of all possible factors affecting the temperature and damage isodose. 



IEEEtran




