
==== Front
ArXiv
ArXiv
arxiv
ArXiv
2331-8422
Cornell University

arXiv:2409.03990v2
2409.03990
2
preprint
Article
Development of Advanced FEM Simulation Technology for Pre-Operative Surgical Planning
Zhao Zhanyue Department of Robotics Engineering, Worcester Polytechnic Institute, Worcester, MA 01605 USA

Jiang Yiwei Department of Robotics Engineering, Worcester Polytechnic Institute, Worcester, MA 01605 USA

Bales Charles Department of Robotics Engineering, Worcester Polytechnic Institute, Worcester, MA 01605 USA

Wang Yang Department of Robotics Engineering, Worcester Polytechnic Institute, Worcester, MA 01605 USA

Fischer Gregory Department of Robotics Engineering, Worcester Polytechnic Institute, Worcester, MA 01605 USA

zzhao4@wpi.edu
9 9 2024
arXiv:2409.03990v2https://creativecommons.org/licenses/by-nc-sa/4.0/ This work is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License, which allows reusers to distribute, remix, adapt, and build upon the material in any medium or format for noncommercial purposes only, and only so long as attribution is given to the creator. If you remix, adapt, or build upon the material, you must license the modified material under identical terms.
nihpp-2409.03990v2.pdf
Intracorporeal needle-based therapeutic ultrasound (NBTU) offers a minimally invasive approach for the thermal ablation of malignant brain tumors, including both primary and metastatic cancers. NBTU utilizes a high-frequency alternating electric field to excite a piezoelectric transducer, generating acoustic waves that cause localized heating and tumor cell ablation, and it provides a more precise ablation by delivering lower acoustic power doses directly to targeted tumors while sparing surrounding healthy tissue. Building on our previous work, this study introduces a database for optimizing pre-operative surgical planning by simulating ablation effects in varied tissue environments and develops an extended simulation model incorporating various tumor types and sizes to evaluate thermal damage under trans-tissue conditions. A comprehensive database is created from these simulations, detailing critical parameters such as CEM43 isodose maps, temperature changes, thermal dose areas, and maximum ablation distances for four directional probes. This database serves as a valuable resource for future studies, aiding in complex trajectory planning and parameter optimization for NBTU procedures. Moreover, a novel probe selection method is proposed to enhance pre-surgical planning, providing a strategic approach to selecting probes that maximize therapeutic efficiency and minimize ablation time. By avoiding unnecessary thermal propagation and optimizing probe angles, this method has the potential to improve patient outcomes and streamline surgical procedures. Overall, the findings of this study contribute significantly to the field of NBTU, offering a robust framework for enhancing treatment precision and efficacy in clinical settings.

Finite Element Modeling
Pre-Operative Surgical Planning
Needle Based Therapeutic Ultrasound
==== Body
pmcI. Introduction

Intracorporeal needle-based therapeutic ultrasound (NBTU) offers a minimally invasive option for the thermal ablation of malignant brain tumors, making it suitable for treating both primary and metastatic cancers [1]. Utilizing a high-frequency alternating electric field (up to 10 MHz) to excite a piezoelectric transducer, NBTU generates acoustic waves that propagate through tissue, causing localized high-temperature heating at the target tumor site. This rapid thermal elevation induces cell death, effectively ablating the tumor. While MR-guided laser interstitial thermal therapy (LITT) has shown experimental success in ablating radio-resistant brain metastases [2], [3], NBTU provides a more precise ablation method by delivering lower doses of acoustic power directly to the tumor or targeted tissue, thereby sparing surrounding healthy tissue with greater efficiency [4]. Moreover, NBTU demonstrates significant potential in treating otherwise inaccessible deep-seated tumors that are near or involve the vascular system [5], [6].

To optimize the energy deposition and element design of NBTU transducers for effective thermal dose delivery during treatment, numerical modeling of the acoustic pressure field generated by the deforming piezoelectric transducer is a critical approach [1]. This modeling is often coupled with simulations of bioheat transfer processes to track the thermal propagation of the applicator over time. Magnetic resonance thermal imaging (MRTI) serves as an experimental validation method for these models, using the proton resonant frequency shift (PRFS) technique to derive quantitative spatial temperature maps from phase differences in the magnetic field. Validation studies using MRTI in gel phantoms have demonstrated the feasibility of these models and their ability to replicate thermal propagation patterns. However, for a more comprehensive evaluation of therapeutic efficacy, thermal damage isodose mapping is preferred. Building on our previous work [1], [7]–[9], this study aims to generate a spectrum of ablation data to support pre-operative surgical planning, providing a robust framework for predicting therapeutic outcomes and enhancing treatment precision.

II. Tissue Damage Estimation

Instead of using a measurement of the amount of energy delivered, thermal dose, which is a measurement of the amount of specific tissue damage caused by heat is widely used. Besides the CEM43 model used in [9], [10], the Arrhenius model is also used in our model. In this model thermal damage is modeled as a first-order rate process in which tissue constituents transform from a native state to a damaged state with a reaction velocity (or rate constant), k. The dependence of k on temperature can be described by Arrhenius Equation 1: (1) k=Ae-EaR(T+273)

Using this assumption the thermal dose, or thermal damage, Ω(τ), is presented as Equation 2 (2) Ω(τ)=lnN(0)N(τ)=A∫0τ e-EaR(T+273)dt

where N(τ) is the concentration of the native tissue at a heating duration of τ, A is the pre-exponential factor (frequency factor), Ea is the activation energy, R is the universal gas constant, and T is the absolute temperature, by adding 273 to make it in degrees of Kelvin. By using the Arrhenius model, tissue damage of both tumor and healthy brain tissue can be evaluated, so that more realistic modeling based on the tumor integrated brain ablation simulation can be achieved.

III. Brain Tumor Integrated Simulation

A. Acoustic Pressure Map in the Trans-Tissue Medium

We used a 3×2mm (a×b) oval shape acting as a tumor inside the brain. In the bioheat physics node additional bioheat material was added and selected tumor oval, then two thermal damage nodes were added under the brain and tumor bioheat material node, the properties parameters used in thermal damage and material properties are shown in Table I. The rest of the configuration was the same as described in [9].

We used a 180° probe and glioma integrated medium as an example, the acoustic pressure pattern can be found in Figure 1, and the results show that the wave propagation was still following the pattern, however, there was a position shift and wave range width increase happened shown in green dashed lines symmetrically, and this was caused by the difference of brain and tumor proprieties. Although there was some misalignment with brain and tumor wave propagation, the pattern was still matching the patterns in our previous work [1], [7], [8], [10], [16]. Based on the results, multiple configurations of medium and parameters simulation were achieved for more complex analysis.

B. Tumor Size Interference

The tumor size had a significant effect on the ablation process, which limited the thermal propagation through the tumor to the surrounding area. Figure 2 shows the temperature propagation through the tumor with different sizes, namely 3×2mm and 8×5mm (a×b) representing the small and large size of glioma, where the properties can be found in Table I, and with the same probe of the same ablation power, and same duration ablation condition. Results show the highest temperature and thermal range was reduced in the small tumor reaching up to 53.081°C versus in the big tumor environment reaching up to 59.303°C. By interfering with the temperature propagation, the thermal damage was influenced as well. Figure 3 shows the CEM43 model damage prediction in different sizes of tumors. The brain tissue and tumor CEM43 map were limited since the temperature was reduced. Because the wave propagation was not changed significantly, no shift or change was observed in the CEM43 of 70 and 100 isodose maps under the trans-tissue condition. However, the area size was changed, which is the same as the reduction of thermal range, in the small tumor the area of CEM43 isodose was reduced compared to the large tumor condition without a trans-tissue scenario.

The Arrhenius damage model was used in the model for further damage estimation, and the result of different size glioma ablation damage estimation can be found in Figure 4. Results show a very bad condition of over-ablated to healthy tissue with the same probe of the same ablation power and same duration ablation condition. Especially in small tumor conditions, the tumor was not ablated yet, while the surrounding healthy tissue was already unexpectedly ablated. In a big tumor, the damage was limited within the tumor range, however the surrounding tissue has started to heat up already, followed by unexpected ablation.

C. Tumor Type Interference

Besides the size of tumor interference, the type of tumor also influences the ablation performance of the different properties between different tumor types. In this study, we compared two common brain tumors, namely glioma and meningioma considering these two types of tumor varies significantly in pre-exponential factor (A). Results of tissue damage are shown in Figure 5, which are also with the same probe of the same ablation power and same duration ablation condition. It is observed that because the meningioma has a larger A factor compared to the glioma, the damage of ablation is more complete, while the glioma is not active enough so the over ablated condition happened in this condition. To avoid a mixed situation, the trans-tissue should be reduced, which indicates that the pre-surgical estimation and planning of ablation should be accurate enough to ablate the target tumor homogeneous medium edge only.

IV. Thermal Ablation Probe Parameter Simulation

We have introduced the tumor integrated ablation simulation, and some major conditions that affected the ablation performance were discussed such as the tumor size and tumor type, which indicated that it is necessary to select the proper probe parameters and there is still room for optimization. In this section, we used the simulation model to achieve the spectrum of ablation performance with parameters combination which provides the database preparation for further studies.

In this study, we chose the four most commonly used directional probes, namely 360°, 180°, 90°, and planar probe with a square tube transducer, and for each probe, we chose 3W, 6W, 10W, 14W and 17W acoustic power under 600s (300s probe on and 300s cooling down) duration time configuration to estimate the trends for our ablation spectrum construction. Following the same analysis procedure, we first achieved the acoustic pressure pattern map, which is shown in Figure 6. Note that all the parameter used for the acoustic medium was based on glioma, and the power was 10W for all the figures shown above. Based on the different power of acoustic pressure pattern maps, we summarized 20 CEM43 model thermal damage isodose maps in glioma, and the results are shown in Figure 7. For different types of probes, the highest temperature reached was proportional to the power used and also varied with probe type. Figure 8 shows the time used to reach the highest temperature for all the simulation conditions with different types of probes and acoustic power. We can observe that the highest temperature reached was related to the probe type under 17W largest acoustic power setup, where the 180° probe reached the highest value of 105.29°C, followed by 90° probe with 87.08°C, 360° probe with 72.78°C, and planar probe with 62.83°C, and reduced power will lower the temperature reached. Moreover, the time used to reach the max temperature was different, where the planar probe used the shortest time of approximately 80s to reach the low increasing period, followed by 90° probe with 90s, 180° probe with 100s, and 360° probe with 150s. All these data are solid preparation for later probe selection protocol construction.

V. Probe Parameters Selection Method for Surgical Pre-Planning

Instead of using the temperature propagation method, the thermal damage model is more important for pre-surgical planning work. In this section, we developed a probe parameters selection method for surgical pre-planning with surgeons and researchers with clinical usage.

Based on the results shown in Figure 7, the pattern of different probes varied a lot with each other, and the selection of directional probe parameters of probe type, power used, and trajectory planning played a vital role in the pre-surgical planning procedure. To address the precise ablation procedure, the following simulation data was accumulated.

Figure 9 and Figure 10 show the CEM43 of 70 and 100 isodose area versus time for all the simulation conditions with different types of probes and acoustic power, which indicates the area of ablation versus time data. This is an area that gradually increases the trend for different probe parameter combinations. Moreover, Figure 11 and Figure 12 show the acoustic power versus CEM43 of 70 and 100 ablated isodose largest area and maximum distance respectively under 600s duration time with four types of probe. This is the most important and valuable database in this work, which acts as a spectrum dictionary and can be used both in the probe selection protocol and further smart ablation studies.

We built a dummy tumor with a 5mm radius as a sample pattern to calculate the approximate ablation time based on all the simulation results we had accumulated. It is observed from Figure 13 that the 360° probe performed the shortest time at approximately 220s compared to other types of probes, namely followed by 180° probe with 300s, 90° probe with 360s, and planar probe with 960s. Figure 14 shows the Ablation angle (θ) generated by four types of directional probes, namely 360° for 360° probe, 60° for 180° probe, 30° for 90° probe, and 5° for the planar probe. Combining the two results above, we can get a summarized characteristic of the four types of probes in this study. The 360° generates a round pattern with the largest pattern range and ablation angle, it is suitable for a round or oval shape tumor and performs the fastest ablation time, however, it is not suitable for the irregular shape of tumors with sharp spikes. On the other hand, the 90° probe generates a flame shape of ablation pattern, which provides the smallest and finest tips of ablation capability, or the highest resolution of ablation, however, it uses the longest time of ablation which is opposite to the 360° probe. The 180° and 90° probe are between the characteristics. A comparison summary of key parameters from four types of probes can be found in Table II. In this study, the time used for ablation is the priority to be considered because a fast procedure improves the patient outcome significantly, so regardless of the power used during the ablation procedure, we suggest choosing a larger degree probe as much as possible as it can be used.

Based on the key parameters compared summary and all the simulation data we have accumulated, a preliminary probe selection protocol for pre-surgical planning towards ablation procedure was developed, which is shown in Figure 15. We first segment the tumor from the scanning images and draw an equivalent circle Ce, this circle is generated by counting the pixels of the tumor portion to get the area value, then using the center of mass as the circle center equivalent O, this equivalent O is the probe insertion target point. The preprocess with tumor is shown in Figure 16. The equivalent radian Re can also be calculated, as well as the internal and external portion compared to the Ce, namely -σ and +σ. Note that the internal portion of -σ can be achieved by decreasing the power during the procedure to reduce the ablation distance. Then some essential parameters are calculated and labeled, namely the maximum distance (Rmax) from the edge with the largest +σ value to equivalent center O following Equation 3, and the ablation angle (θ). Note that for compatibility consideration, we choose the minimum ablation angle for this step, namely θmin. The next step is the acoustic power selection, which is based on the Rmax value and looks for compatible acoustic power from all four types of probes from the power versus CEM43 maximum isodose distance data set (Figure 12). The last step is to finalize the probe, which is based on the θmin and the CEM43 isodose pattern propagation angle range (Figure 14). In this step, we will choose the probes with an ablation range angle that is smaller than θmin so that no over ablation issue will happen. Note that if multiple probes are still under consideration after all the protocols, choose the largest ablation angle range probe for the final decision. (3) Rmax=Re+(+σ)

Figure 17 shows 3 samples of tumor after labeled with Re,±σ portions, and θmin. (a) the tumor is almost a circle with very small ±σ portions, so this tumor is suitable with a 360° probe for rapid ablation. (b) the tumor is a H2 molecule shape, although labeled with essential parameters, and the θmin is larger than 60°, which means it is suitable with 180° probe, however, this pattern is more suitable for multiple entry points insertion with multiple ablation duration. This is not in our study and requires further optimization study. And lastly (c) tumor is an irregular shape with two θ values namely approximately 45° and 10°. Considering the 10° of +σ portion is very small and the edge is almost located on the edge without a large spike compared to the 45° portion, we can choose the 90° probe for rapid ablation or the planar probe for the high resolution ablation procedure.

VI. Conclusion and Disscussion

This work presents the database preparation for future studies based on the extended simulation results. We first developed an extended simulation integrated with different tumors to evaluate the effect of thermal damage in trans-tissue conditions. Then we accumulated the simulation database of four directional probes, including the CEM43 isodose map, temperature change, thermal dose area, and maximum distance. Finally, an idea of the probe selection method was presented for pre-surgical planning.

In the tumor integrated simulation, we involved the Arrhenius equation for tissue damage estimation, and results showed a significant effect from the tumor size and type with the trans-tissue scenario, so we suggested avoiding thermal propagation and damage model with trans-tissue modeling, as well as in the real clinical surgery. Also, the change of parameters during the ablation process should also be considered, especially when dealing with multiple tissues in a complex environment.

Then the database of four directional probes, including the CEM43 isodose map, temperature change, thermal dose area, and maximum distance based on the simulation model built up the dictionary of thermal ablation therapy, this work will benefit the future complex trajectory planning. Furthermore, this simulation method can also build up an easy-to-use platform for researchers to develop various parameters in combination with NBTU thermal ablation study.

Finally, the idea of the directional probe selection method gave us inspiration for how to use the database developed above, and this method will also benefit the pre-surgical planning study and eventually reduce the study curve for surgeons and researchers. The idea of choosing a large angle probe will significantly improve the outcomes of patients by reducing the ablation procedure time.

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

Fig. 1. Acoustic pressure pattern of a 180° probe and glioma integrated medium.

Fig. 2. Thermal propagation maps in different sizes of tumors. The highest temperature reached is influenced by the size of the tumor.

Fig. 3. CEM43 isodose map in different sizes of tumors. The area of CEM43 70 and 100 is influenced by the size of the tumor.

Fig. 4. Tissue damage is estimated in different sizes of tumors. The small tumor is not ablated yet, the surrounding tissue already died, while for the big tumor, the damage was limited within the tumor range, however the surrounding tissue has started to heat up already, followed by unexpected ablation.

Fig. 5. Tissue damage is estimated in different types of tumors, namely glioma and meningioma.

Fig. 6. Acoustic pressure pattern from the 4 most common directional probes used in the study.

Fig. 7. CEM43 isodose maps of 360°, 180°, 90°, and planar probe under 3W, 6W, 10W, 14W, and 17W acoustic power respectively.

Fig. 8. The highest temperature versus time for all the simulation conditions with different types of probes and acoustic power.

Fig. 9. The CEM43 of 70 isodose areas versus time for all the simulation conditions with different types of probes and acoustic power.

Fig. 10. The CEM43 of 100 isodose area versus time for all the simulation conditions with different types of probes and acoustic power.

Fig. 11. Acoustic power versus CEM43 of 70 and 100 largest ablated isodose area under 600s duration time with four types of probe.

Fig. 12. Acoustic power versus CEM43 of 70 and 100 maximum ablated isodose distance under 600s duration time with four types of probe.

Fig. 13. Ablation complete-time was used from four types of directional probes, under a 5mm radian round sample tumor pattern.

Fig. 14. Ablation angle (θ) generated by four types of directional probes, namely 360° for 360° probe, 60° for 180° probe, 30° for 90° probe, and 5° for the planar probe.

Fig. 15. Probe selection protocol for pre-surgical planning.

Fig. 16. Segment the tumor scan images, draw the equivalent circle, and label the circle center.

Fig. 17. 3 tumor samples after labeled with Rmax,±σ portions, and θmin

TABLE I Acoustic medium properties of healthy brain tissue, glioma, and meningioma. GM - gray matter. Data imported from [11]–[15].

Acoustic Medium Properties	Units	Brain (GM)	Glioma	Meningioma	
Heat Capacity (C)	J / kg/° C	3680	3680	3680	
Density (ρ)	kg/m 3	1035.5	1046	1058	
Thermal Conductivity (K)	W/m/° C	0.565	0.565	0.565	
Blood Capacity (Ct,b)	J / kg/° C	3630	3630	3630	
Blood Density (ρb)	kg/m 3	1050	1050	1050	
Blood Perfusion (ωb)	kg/m3 / s	0.013289	0.0005	0.0005	
Speed of Sound (c)	m/s	1460	1660	1590	
Metabolic Heat (Qm)	W/m 3	16229	25000	25000	
Frequency Factor (A)	l/s	1.66e91	3.153e47	1.737e91	
Activation Energy (Ea)	J/mole	5.68e5	3.149e5	3.149e5	

TABLE II Key parameters of four types of probes compare summary.

Probe	Range	Time	Resolution	
360°	Large	Fast	Low	
180°	Moderate	Mid-Fast	Moderate	
90°	Mid-Small	Moderate	Mid-high	
Planar	Small	Slow	High
==== Refs
References

[1] Gandomi K. , Zhao Z. , Tarasek M. , Fiveland E. , Bhushan C. , Ghoshal G. , Neubauer P. , Williams E. , Liu M. , Carvalho P. , “3d thermo-acoustic modeling of a piezoelectric transducer for directed interstitial ultrasound ablation.”
[2] Rao M. S. , Hargreaves E. L. , Khan A. J. , Haffty B. G. , and Danish S. F. , “Magnetic resonance-guided laser ablation improves local control for postradiosurgery recurrence and/or radiation necrosis,” Neurosurgery, vol. 74 , no. 6 , pp. 658–667, 2014.24584138
[3] Carpentier A. , McNichols R. J. , Stafford R. J. , Guichard J.-P. , Reizine D. , Delaloge S. , Vicaut E. , Payen D. , Gowda A. , and George B. , “Laser thermal therapy: Real-time mri-guided and computer-controlled procedures for metastatic brain tumors,” Lasers in surgery and medicine, vol. 43 , no. 10 , pp. 943–950, 2011.22109661
[4] Ghoshal G. , Salgaonkar V. , Wooton J. , Williams E. , Neubauer P. , Frith L. , Komadina B. , Diederich C. , and Burdette E. C. , “Ex-vivo and simulation comparison of multi-angular ablation patterns using catheter-based ultrasound transducers,” in Energy-based Treatment of Tissue and Assessment VII, vol. 8584 . SPIE, 2013, pp. 293–303.
[5] Missios S. , Bekelis K. , and Barnett G. H. , “Renaissance of laser interstitial thermal ablation,” Neurosurgical focus, vol. 38 , no. 3 , p. E13, 2015.
[6] Christian E. , Yu C. , and Apuzzo M. L. , “Focused ultrasound: relevant history and prospects for the addition of mechanical energy to the neurosurgical armamentarium,” World Neurosurgery, vol. 82 , no. 3–4 , pp. 354–365, 2014.24952224
[7] Gandomi K. , Carvalho P. , Zhao Z. , Nycz C. , Burdette E. , and Fischer G. , “Thermo-acoustic simulation of a piezoelectric transducer for interstitial thermal ablation with mrti based validation,” in Comsol Conference, 2019.
[8] Gandomi K. Y. , Carvalho P. A. , Tarasek M. , Fiveland E. W. , Bhushan C. , Williams E. , Neubauer P. , Zhao Z. , Pilitsis J. , Yeo D. , “Modeling of interstitial ultrasound ablation for continuous applicator rotation with mr validation,” IEEE Transactions on Biomedical Engineering, vol. 68 , no. 6 , pp. 1838–1846, 2020.
[9] Zhao Z. , Szewczyk B. , Tarasek M. , Bales C. , Wang Y. , Liu M. , Jiang Y. , Bhushan C. , Fiveland E. , Campwala Z. , “Deep brain ultrasound ablation thermal dose modeling with in vivo experimental validation,” arXiv preprint arXiv:2409.02395, 2024.
[10] Szewczyk B. , Tarasek M. , Campwala Z. , Trowbridge R. , Zhao Z. , Johansen P. M. , Olmsted Z. , Bhushan C. , Fiveland E. , Ghoshal G. , “What happens to brain outside the thermal ablation zones? an assessment of needle-based therapeutic ultrasound in survival swine,” International Journal of Hyperthermia, vol. 39 , no. 1 , pp. 1283–1293, 2022.36162814
[11] Sadeghi-Goughari M. , Mojra A. , and Sadeghi S. , “Parameter estimation of brain tumors using intraoperative thermal imaging based on artificial tactile sensing in conjunction with artificial neural network,” Journal of Physics D: Applied Physics, vol. 49 , no. 7 , p. 075404, 2016.
[12] de Oliveira M. M. , Wen P. , and Ahfock T. , “Heat transfer due to electroconvulsive therapy: Influence of anisotropic thermal and electrical skull conductivity,” Computer methods and programs in biomedicine, vol. 133 , pp. 71–81, 2016.27393801
[13] Vaupel P. , Kallinowski F. , and Okunieff P. , “Blood flow, oxygen and nutrient supply, and metabolic microenvironment of human tumors: a review,” Cancer research, vol. 49 , no. 23 , pp. 6449–6465, 1989.2684393
[14] Bousselham A. , Bouattane O. , Youssfi M. , and Raihani A. , “Brain tumor temperature effect extraction from mri imaging using bioheat equation,” Procedia Computer Science, vol. 127 , pp. 336–343, 2018.
[15] Tanaka K. , Ito K. , and Wagai T. , “The localization of brain tumors by ultrasonic techniques: A clinical review of 111 cases,” Journal of Neurosurgery, vol. 23 , no. 2 , pp. 135–147, 1965.5842309
[16] Campwala Z. , Szewczyk B. , Maietta T. , Trowbridge R. , Tarasek M. , Bhushan C. , Fiveland E. , Ghoshal G. , Heffter T. , Gandomi K. , “Predicting ablation zones with multislice volumetric 2-d magnetic resonance thermal imaging,” International Journal of Hyperthermia, vol. 38 , no. 1 , pp. 907–915, 2021.34148489
