Natural bond orbital analysis of dication magnesium complexes [Mg(H₂O)₆]²⁺ and [[Mg(H₂O)₆](H₂O)_n]²⁺; n=1-4

1.   Introduction

With the rising incidence of cancer, predicted to triple by 2030 according to the World Health Organization (WHO)^[1], breast cancer remains a major concern, accounting for 284, 200 new cases and 44, 130 deaths in the United States in 2021^[2]. The limitations of surgical interventions, particularly for early-stage tumors, have fueled demand for less invasive treatments, with the focused ultrasound (FUS) emerging as a promising alternative. FUS precisely targets and destroys cancer cells using ultrasound waves, offering noninvasiveness, cost-effectiveness, and enhanced drug delivery capabilities ^[3,4,5]. Recent studies, such as those by Deckers et al.^[6], Peek et al.^[7], and Feril et al. ^[8], underscore FUS's potential in breast cancer treatment, yet there remains a significant gap in comprehensive studies simulating its complex multi-physics phenomena on realistic models. Our research addresses this by developing a deep-learning model to assess FUS's effectiveness, building on recent analytical and numerical insights.

In recent years, significant advancements have been made in simulating the effects of FUS on various tissues. Yoon and colleagues ^[9] developed a multi-resolution simulation using the finite-difference time-domain formulation for acoustic wave transmission, while Gupta et al.^[10] and Rezaeian et al.^[11] investigated how tissues absorb and dissipate heat during FUS, and the delivery of thermosensitive liposomal doxorubicin using FUS, respectively. Montienthong et al.^[3] conducted a simulation focusing on FUS ablation in breast cancer treatment using a two-dimensional model, while Gupta et al.^[12] and Kaczmarek et al.^[13] explored bio-heat transfer models, including the thermal wave model and Pennes' bio-heat equation. Mohammadpour et al.^[14] further studied heat transfer and hemodynamics in FUS ablation of the porous liver. However, the literature on FUS is constrained by the lack of comprehensive 3D modeling and anatomically accurate breast phantoms. To address this, machine learning and deep neural networks, particularly physics-informed neural networks (PINNs), have been increasingly leveraged in scientific computing, providing an alternative to traditional numerical methods and enabling the simulation of realistic clinical procedures and detailed examinations of FUS ablation processes ^{[15,16,17,18,19]}.

PINNs have revolutionized the approach to solving nonlinear partial differential equations (PDEs) by leveraging deep learning to approximate solutions while adhering to physical laws, significantly reducing computational demands compared to traditional methods like the finite difference method (FDM) and finite element method (FEM) ^{[20,21,22,23,24]}. This study utilizes a 3D model based on an anatomically realistic breast phantom (ARBP) generated from magnetic resonance imaging (MRI), accurately reflecting the breast's complex anatomy and enhancing the precision and clinical relevance of FUS simulations ^[25]. The research explores realistic FUS ablation processes with multiple lesion points, assessing their cumulative effect on necrotic tissue, demonstrating the deep-learning platform's ability to mirror actual clinical treatments ^[26]. The paper presents a comprehensive method for predicting FUS therapy outcomes using an ARBP, detailing a multi-stage simulation process that includes generating an acoustic wave, computing acoustic intensity distribution, and employing the absorbed acoustic energy in a bio-heat transfer model to estimate the generated heat and resulting temperature increase in the tumor area. This allows for determining the percentage of necrotic tissue. The paper is organized as follows: Section 2 details the simulation process, Section 3 verifies the accuracy of the proposed simulation, Section 4 discusses deep learning results concerning temperature elevation and necrotic tissue percentage, including a validation against FEM outcomes, and Section 5 concludes the study with key findings.

2.   Materials and methods

This section provides a concise overview of the FUS ablation process. It proceeds to elaborate on the governing equations and mathematical models involved. The following subsections offer vital details concerning the deep learning modeling procedure at various stages.

2.1.   FUS ablation mechanism

FUS ablation is a noninvasive procedure where an ultrasound transducer generates focused ultrasonic beams that penetrate soft tissue and specifically target tumors ^[27]. The focused beams concentrate high intensities within a compact volume, typically around 1 mm in diameter and 10 mm in length, with sound intensity in the focal region ranging from 100 to 10, 000 Wcm^-2 and peak pressures reaching up to 70 MPa (compression) and 20 MPa (rarefaction). These high levels of acoustic intensity and pressure are primarily due to nonlinear effects, which significantly enhance energy concentration at the focal point. The two primary mechanisms involved in FUS ablation are the thermal effect and mechanical effect ^[28]. The thermal effect involves heat generation from acoustic energy absorption, leading to a rapid temperature increase in the local tissue. When the temperature exceeds 60 ℃ for just 1 second, it typically causes rapid and irreversible cell death through coagulation necrosis, which is the primary mechanism for eliminating tumor cells during FUS treatment ^[28].

The mechanical effects induced by FUS, including cavitation, micro-streaming, and radiation force, specifically occur when high-intensity acoustic pulses are employed. However, due to the thermal process being better understood and more controllable, it is preferred for tissue ablation ^[28]. FUS is rapidly gaining recognition in the medical field as a noninvasive ablation technique across various applications. Treatments are typically completed in a single session, with the patient being fully conscious, lightly sedated, or under mild general anesthesia ^[29]. To enhance comprehension, the schematic depiction of the FUS ablation process for a breast tumor is illustrated in Figure 1.

Figure 1.  A schematic of the FUS ablation process for the treatment of a breast tumor.

DownLoad: Full-Size Img PowerPoint

2.2.   Governing equations

In this section, we examine the essential equations and steps necessary to comprehensively simulate the process of FUS ablation. Initially, we solve the equations governing the behavior of acoustic waves as they propagate within the breast tissue. Subsequently, we simulate the generation of heat resulting from the absorption of these acoustic waves and the subsequent transfer of heat within the biological tissue of the breast. Ultimately, we make predictions regarding the extent of tumor damage and the occurrence of necrosis. Presented below are the governing equations that need to be addressed in a specific order during the solution process.

2.2.1.   Acoustic wave generation and propagation

To represent the fixed acoustic field, we address the inhomogeneous Helmholtz equation through our modeling approach. The specific form of this equation can be found in reference ^[30]:

In the equation, $\rho$ represents density, $p$ represents total pressure, $k$ represents the wave number, $q_d$ represents the source of the dipole domain, and $Q_m$ represents the source of the monopole domain. As acoustic waves travel through breast tissues, they experience a reduction in intensity known as attenuation. Due to the involvement of different tissue types, a customized attenuation model is employed for each tissue, utilizing specific attenuation coefficients ($\alpha_a$). This inclusion of attenuation results in the wave number becoming a complex value, as depicted in the following equation ^[30]:

In the equation, $\omega$ represents the angular frequency and $c$ represents the speed of sound.

2.2.2.   Bio-heat transfer in breast tissues

In this context, the simulation of temperature distribution in biological tissues is achieved by employing Penne's bio-heat transfer equation, as presented in reference ^[31]:

The mathematical equation involves multiple variables corresponding to the physical properties of both tissue and blood, where $\rho$ represents tissue density, $c_p$ is the specific heat at constant pressure, $k_h$ is the heat conduction coefficient, $\omega_b$ indicates the local blood perfusion rate, $\rho_b$ is the blood density, $c_{p, b}$ is the blood's specific heat, $T_b$ refers to the arterial blood temperature, $Q_{met}$ signifies the metabolic heat generation rate, and $Q_S$ denotes the heat source rate corresponding to the power absorbed by the acoustic waves. Once the acoustic field is resolved, the acoustic pressure and intensity fields are determined, allowing the calculation of the $Q_S$ value required for thermal simulation according to references ^[30,32].

The equation mentioned $\alpha_a$ represents the attenuation coefficient, $I$ represents the magnitude of acoustic intensity, $p$ represents the acoustic pressure, and $v$ represents the vector of acoustic particle velocity. Additionally, the magnitude of sound intensity refers to the power conveyed by the acoustic waves per unit area in a direction perpendicular to that area. It is defined as stated in reference ^[33]:

The equation, $p_A$ represents the amplitude of pressure oscillation, while $\rho$ and $c$ represent the density and speed of sound, respectively. The pressure distribution should be considered along the entire wavefront to account for the maximum acoustic pressure. This is crucial for accurately determining the intensity of the ultrasound wave. Moreover, acoustic intensity is defined as the power per unit area carried by the wave. The intensity of an ultrasound wave is directly proportional to the square of the acoustic pressure. Therefore, using the peak pressure value in the calculation provides an accurate estimation of the maximum possible intensity, which is essential for effective and safe FUS treatment.

2.2.3.   Prediction of the fraction of necrotic tissue

The evaluation of tissue injury resulting from the hyperthermia process can be conducted using the Arrhenius kinetic model, as described in references ^[34,35]:

In the equation mentioned above, $n$ denotes the polynomial order within the equation. $R$ represents the universal gas constant, while $T$ corresponds to the tissue temperature. $A$ and $\Delta E$ are frequency factors and activation energy, respectively, which are specific to the particular tissue type and have been previously determined for various tissue types. For breast tissue, these parameters were calculated as $A$ = 1.18 × 10⁴⁴ s^-1 and $\Delta E$ = 3.02 × 10⁵ Jmol^-1. Once the degree of tissue injury ($\alpha$) is determined, the fraction of necrotic tissue ($\theta_d$) can be expressed accordingly, as detailed in reference ^[35]:

2.3.   Geometry of the model

All geometry processing in this article utilizes the commercial software "MATLAB" ^[36]. The ARBP phantom is placed inside a water tank, with an acoustic transducer submerged at the bottom to generate an ultrasound wave that propagates through the water and breast tissues. The transducer, shaped like a bowl with a focal length of 51.74 mm, is aligned with the tumor's center to focus acoustic pressure, leading to rapid temperature increase and tissue ablation, as depicted in Figure 2(a). A cross-sectional view of the model, shown in Figure 2(b), highlights the transducer's focal length and aperture of 51.74 and 64 mm, respectively. This section details the use of 3D grid-based numerical breast models with accurate anatomical details in deep learning and FEM simulations, utilizing a model from the University of Wisconsin Cross-Disciplinary Electromagnetics Laboratory (UWCEM) repository ^[25]. These models, derived from T1-weighted MRIs of patients in a face-down position, are designed to accurately represent breast anatomy for cancer detection. The selected model, Breast ID 070604PA1 (Class 2), was chosen based on the American College of Radiology's classification system, which categorizes breasts by radiological density, with Class 2 representing some fibro-connective/glandular tissue ^[37].

Figure 2.  (a) 3D visualization of the entire geometry considered; (b) Middle cross-section view of the geometry.

DownLoad: Full-Size Img PowerPoint

By utilizing the previously mentioned data and a MATLAB ^[36] function connected to the breast phantom's tissue types, a 3D detailed surface representation is generated within spatial coordinates. These surfaces, representing different tissue types, are processed using SOLIDWORKS software to convert them into surfaces and volumes, resulting in a comprehensive 3D model of the breast phantom ^[38]. The process involves connecting relevant points from each tissue type to form surfaces, which are then interconnected to create a closed volume that accurately replicates the breast phantom's tissue types. To simplify this, four specific tissue types—fibro-connective and glandular tissue (FCG), transitional tissue, fatty tissue, and muscle—are isolated from the database. Figure 3 illustrates the steps involved in generating a 3D breast phantom, from the detailed surface to the resulting volume corresponding to fatty tissue. In our simulations, the tumor is represented as an ellipsoid object with a 10 mm z-diameter and 5 mm x and y diameters, classified as T1c, indicating its small size and early stage of cancer development ^[39]. The tumor, located within the fatty tissue, along with the other tissue types, is collectively represented as the ARBP in Figure 4.

Figure 3.  Three distinct processes were used to build a 3D breast phantom: (a) obtaining a detailed surface from the database; (b) constructing surfaces from the detailed surface; and (c) producing the volume from knitted surfaces.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Considered 3D geometry of tissue types separately and collectively: (a) FCG tissue; (b) transitional tissue; (c) fatty tissue, muscle, and tumor; and (d) all tissues within each other as the ARBP.

DownLoad: Full-Size Img PowerPoint

2.4.   Simulation procedure

This section provides a comprehensive explanation of the numerical modeling process for FUS, detailing the necessary data, material properties, and initial and boundary conditions for each physics involved. The process begins with simulating the generation of ultrasonic waves within the specified domain, followed by heat generation and transfer within the breast phantom. The initial step focuses on pressure acoustics in frequency-domain physics, where an input displacement amplitude ($d_0$) of 3.8 nm is introduced at the transducer aperture boundary. According to Diaz et al. ^[40], solutions remain within the linear regime when the input pressure is below $p_0$ = 0.35 MPa, but deviate into a nonlinear regime as pressure increases, especially in the focal region. The solution assumes linear propagation of acoustic waves, neglecting nonlinear effects and shear waves, with an input pressure of around 0.1 MPa and an acoustic frequency of 1.1 MHz.

Following the simulation of acoustic wave propagation within the breast phantom, the next step involves incorporating time-dependent bio-heat transfer physics. The acoustic intensity field is determined based on the acoustic pressure field, and the heat source ($Q_S$) is calculated using Eq (2.4). In our simulations, the acoustic waves are present for 1 second, with the ablation process lasting 5 seconds. The heat source generated by the ultrasonic waves remains active for only 1 second, after which the simulation continues without it for an additional 4 seconds to simulate the cooling process and determines the extent of necrotic tissue resulting from the 5-second exposure to elevated temperatures. Thermal damage is calculated using the Arrhenius kinetics formula (Eq (2.5)), with considerations for heat transfer through free convection, radiation, and the impact of blood perfusion on cooling the breast and tumor area, as outlined in Eq (2.3). The properties of blood, blood perfusion rate, and metabolic heat source for each tissue type are provided in Table 1, and Figure 5 illustrates the boundary conditions for each physics.

Table 1.  Thermophysical and acoustic properties of different materials ^{[1,41,42,43,44,45,46,47,48]}.

Materials	$\rho$	$C_p$	$k$	$c$	$\alpha_a$	$Q_{met}$	$\omega_b$
	(kgm-3)	(Jkg-1K-1)	(Wm-1K-1)	(ms-1)	(dBm-1MHz-1)	(Wm-3)	(s-1)
FCG tissue	1050	3770	0.48	1470	87	700	0.0067
Transitional tissue	990	3270	0.345	1463.5	75	700	0.0067
Fatty tissue	930	2770	0.21	1457	48	700	0.0067
Tumor	1050	3852	0.54	1509	57	5790	0.005
Muscle	1100	3800	0.48	1588.4	109	5790	0.005
Water	997	4184	0.598	1480	0.22	-	-
Blood	1050	3617	0.52	1540	20	-	-

---

 | Show Table

DownLoad: CSV

Figure 5.  Boundary conditions considered in different physics, (a) pressure acoustics and (b) bio-heat transfer.

DownLoad: Full-Size Img PowerPoint

As depicted in Figure 5, the boundary conditions for frequency-domain pressure acoustics involve considering a sound hard boundary for the water tank and transducer walls. Given that the chest wall and other boundaries within the breast tissue are far from the acoustic source (the transducer surface), they are also assumed to be sound hard boundaries. While the transducer surface is modeled as a hard boundary, this does not imply complete rigidity; rather, it means that the normal components of particle velocity and acceleration are zero at the boundary, except at the transducer surface, which is designed to vibrate and emit ultrasonic waves. Although perfect matching layers (PML) could be advantageous in absorbing outgoing waves and reducing reflections, they increase boundary condition complexity and computational demands, especially in heterogeneous media like breast tissue. Therefore, a sound-hard boundary is preferred in the modeling process. To apply the input pressure signal, the transducer surface is assigned a normal displacement boundary condition, with an input displacement amplitude of 3.8 nm. For heat transfer modeling, bio-heat transfer physics incorporates both heat flux and radiation from the skin to the surrounding water as boundary conditions. The temperature across internal boundaries is assumed to be the normal human body temperature of 37℃, set as the temperature boundary condition for accuracy, as detailed in the corresponding equations ^[35].

The equations in this study involve various variables: $\overrightarrow{n}$ represents the surface normal vector, $\rho$ stands for density, $p$ represents total pressure, $q_d$ represents the dipole domain source, $d_0$ is the displacement amplitude of the input pressure signal, $q$ represents the heat transfer vector, $h$ is the heat convection coefficient, $T_\infty$ represents the ambient temperature, $T$ represents the temperature, $\epsilon$ is the surface emissivity, and $\sigma_{SB}$ represents the Stefan-Boltzmann constant. The heat convection coefficient between the breast skin and surrounding water is set to 53.5 Wm^-2K^-1, based on the study by Boutelier et al. ^[49], with an ambient air temperature of 293.15 K and a skin surface emissivity of 0.98 ^[50]. Detailed material properties used in the study are provided in Table 1.

The provided table contains various properties such as density ($\rho$), specific heat under constant pressure ($C_p$), heat conductivity ($k$), speed of sound ($c$), acoustic attenuation coefficient ($\alpha_a$), metabolic heat generation rate ($Q_{met}$), and blood perfusion rate ($\omega_b$). To assess the simulation outcomes, three cross-sectional planes were taken into account, intersecting at the center of the tumor. These planes are depicted in Figure 6 and consist of the x-y, x-z, and y-z surfaces.

Figure 6.  Considered cut planes that pass through the center of the tumor.

DownLoad: Full-Size Img PowerPoint

2.5.   PINN

Training a deep learning algorithm to accurately identify a nonlinear mapping between high-dimensional input and output data pairs is challenging. Researchers have found that incorporating empirically validated rules as regularization factors can guide time-dependent dynamics within a predefined system, limiting the solution space to a permissible domain and enhancing the algorithm's ability to achieve optimized solutions and generalize well, even with limited training samples ^{[23,24,51,52,53]}. Recent progress in physics-based deep learning techniques has demonstrated that neural networks, as illustrated in Figure 7(a), can effectively handle PDEs by integrating them into the loss function along with initial conditions, boundary conditions, and other constraints, enabling the network to predict the behavior of a computational domain governed by specific physics laws.

Figure 7.  (a) Diagram of PINN architecture for solving PDE. The input layer of the neural network with spatial inputs and temporal. The output layer computes the solution of the given PDE and embeds it into the loss function. PDE can have constraints such as initial conditions and boundary conditions. By minimizing the loss function, PINN can generate the solution with respect to the previous conditions; (b) diagram of PINN architecture for solving the 3D Helmholtz equation and time-dependent bio-heat transfer equation (tBHTE) to predict the temperature rise. The upper neural network with three spatial inputs $x$, $y$, and $z$, and the pressure result $p$ as the output value. The lower neural network takes three previous spatial inputs as well as a temporal input, $t$, and produces the temperature result, $T$, as the output value. Other sections are PDE, initial conditions (IC) and boundary conditions (BC). The last section represents the loss functions embedded with PDE, BC, and IC along with the training output $T$. By minimizing the loss function, PINN can generate the heat map with respect to the previous conditions.

DownLoad: Full-Size Img PowerPoint

PINNs integrate physical laws into neural networks by embedding governing PDEs directly into the loss function, which consists of two main components: data loss and physics loss. The data loss measures the difference between predicted outputs and actual data, ensuring accurate approximation, while the physics loss penalizes deviations from the PDEs, expressed through the network's outputs and their derivatives using automatic differentiation. This approach balances fitting observed data with adhering to physical laws, enhancing generalization and ensuring physically consistent predictions. In the deep learning configuration, waves from concave ultrasound transducers are modeled by solving the Helmholtz equation and used as input for the next PINN simulator. The predicted output then serves as an external heat source in the tBHTE, which is fed into a second PINN structure that predicts heat conduction and generates a temperature map within the tissue domain.

The hyperparameters for training the PINN model were meticulously chosen to optimize performance, with a learning rate of 0.001 selected based on preliminary experiments to ensure stable and efficient training—lower rates slowed convergence, while higher rates caused instability. A batch size of 32 balanced memory usage and training speed, and training was conducted over 5000 epochs, determined by monitoring the validation loss until it plateaued, indicating model convergence. The Adam optimizer was used for its efficiency and adaptive learning rate capabilities. For the plot illustrating wave propagation and the heat map, the solution of PDEs using PINN involved setting up components like temporal and spatial domains, the frequency-domain Helmholtz equation, the tBHTE, boundary conditions, initial conditions, and other hyperparameters. The spatial and temporal domains were set at 1 mm in each direction and 0.1 s, respectively, with the model trained to ensure that given an input $(x, y, z, t)$, the outputs $p(x, y, z)$ and $T(x, y, z, t)$ satisfy the wave and heat equations along with the boundary and initial conditions.

The loss function for the PINN model is defined by two components: $MSE_f = \frac{1}{N}\sum_{i = 1}^{N_f}|f(x_f^i, y_f^i, z_f^i, t_f^i)|^2$, representing the mean squared error over the interior domain points $[x_f^i, y_f^i, z_f^i, t_f^i]_{i = 1}^{N_f}$, and $MSE_u = \frac{1}{N}\sum_{i = 1}^{N_u}|f(x_u^i, y_u^i, z_u^i, t_u^i)-u^i|^2$, which accounts for the BCs and ICs at the collocation points $[x_u^i, y_u^i, z_u^i, t_u^i]_{i = 1}^{N_u}$. The neural network model employs 9 hidden layers with 20 neurons each, a learning rate of 0.001, and the rectified linear unit (ReLU) activation function. Hyperparameter tuning is crucial, as different settings may require adjustments to capture the underlying physics accurately. Increasing the number of layers or neurons can improve the model's ability to learn complex patterns but may risk overfitting if not properly regularized, while a simpler model might fail to capture essential features. Therefore, careful tuning of hyperparameters like the number of layers, neurons, learning rate, and activation functions is essential for optimizing model performance and achieving a balance between complexity and accuracy.

Table 1 presents the essential attributes for simulating the temperature increase caused by a concave ultrasound transducer, including the parameters related to the tBHTE and the external heat source $Q_S$. After implementing the wave equation and constraints for creating a concave transducer in the neural network (Eq (2.1)), the resulting value is used as the external heat source in the simulation. This allows the PINN to model heat conduction while considering all the relevant features. Subsequently, the tBHTE is incorporated into the neural network as the PDE to solve the problem outlined in Eq (2.3).

Subsequently, our model incorporates the embedded equation within the neural network, which has 9 hidden layers with 20 neurons in each layer, to forecast the output of the associated tBHTE using spatial and temporal inputs. The training output is responsible for simulating thermal conduction within a real-life bio-tissue domain. By employing a deep learning algorithm to approximate the solution of the PDE, the neural network minimizes the residual. The depicted neural network in Figure 7(b) represents the architecture of the PINN utilized in the proposed methodology.

In contrast to traditional machine learning models, the PINN approach does not require separate validation and testing datasets. The validation process is inherent in the model's structure, as it is guided by the governing PDEs, which are embedded in the loss function. These PDEs ensure that the model's predictions are physically consistent and accurate during training. For testing purposes, we compare the PINN model's results with those from the FEM. This comparison allows us to quantitatively assess the model's accuracy in predicting temperature distribution, acoustic pressure, and necrotic tissue formation under various clinical scenarios. The alignment between the PINN and FEM results confirms the robustness and reliability of the PINN approach.

2.6.   Finite element modeling

The FEM simulation utilizes quadratic Lagrange discretization for both the pressure acoustics and bio-heat transfer phenomena. The FEM simulation was performed using Comsol Multi-physics software on a system equipped with an Intel Core i7-1065G7 Processor running at 1.5 GHz and 16 GB DDR4 RAM. Each simulation took approximately 115 hours to complete. The majority of the computation time is dedicated to simulating the acoustic wave due to the need for very fine meshes in this particular physics simulation. The subsequent subsection will provide a detailed explanation of the required mesh grids for each step mentioned. The various boundary conditions considered for each of the mentioned physics can be observed in Figure 5.

2.6.1.   Mesh study

The mesh size is critical in accurately simulating acoustic waves, particularly in FEM simulations of "pressure acoustics, " which require significantly finer meshes compared to "bio-heat transfer" simulations. To achieve precise solutions for acoustic equations, a highly refined mesh set with a maximum mesh size of $\lambda/4$ (where $\lambda$ represents the wavelength) is employed for the ARBP and the surrounding water. For the tumor tissue, an even finer mesh with a maximum size of $\lambda/8$ is used to capture the sharp pressure gradient in the focal region accurately. Additionally, finer meshes are applied near the focal point, while larger mesh sizes are considered for regions farther from the focal point, as they have a lesser impact on the accuracy of the results.

In this research, multiple computational meshes were generated to assess grid independence by computing the maximum sound intensity magnitude at the focal point. The ideal number of grids was determined by balancing computational cost and simulation accuracy, with the most refined grid selected if finer meshes did not significantly alter the results, as shown in Figure 8(a) for "pressure acoustics" physics. For solving bio-heat transfer equations, significantly larger meshes were used to achieve high precision while minimizing simulation time. Similar mesh samples were analyzed to measure the maximum temperature at the focal point during the ablation process, with the grid independence results for "bio-heat transfer" physics presented in Figure 8(b).

Figure 8.  Grid independency test for each physics.

DownLoad: Full-Size Img PowerPoint

Both Figure 8(a), (b) show that as the number of grids increases, the maximum intensity magnitude and temperature at the focal point stabilize. Table 2 lists numerical values and differences across simulated scenarios to identify the optimal mesh for "pressure acoustics" and "bio-heat transfer" physics. Relative errors, which can be positive (indicating higher measured values) or negative (indicating lower values), reveal that the maximum intensity magnitude for cases 1.3 and 1.4 is nearly identical, making case 1.3 the optimal mesh for "pressure acoustics" to reduce computational costs. Similarly, cases 2.4 and 2.5 show comparable maximum temperatures, with case 2.4 chosen as the optimal mesh for "bio-heat transfer." Figure 9(a), (c) display these optimal mesh designs, while sliced views in Figure 9(b), (d) provide a clearer visualization of the individual mesh grids for each physics.

Table 2.  Quantitative amounts and relative errors of different considered mesh cases for each of the physics.

Pressure acoustic physic
Case	Grid numbers	Maximum sound intensity magnitude (Wm-2)	% Relative error in comparison with the optimum case
1.1	178438795	5276891	-7.471 %
1.2	359274873	5560187	-2.504 %
1.3	447531449	5702988	-
1.4	519592973	5702988	0.116 %
Bio-heat transfer physic
Case	Grid numbers	Maximum temperature generated at focal point (℃)	% Relative error in comparison with the optimum case
2.1	1421843	62.90969	-3.788 %
2.2	1798417	64.14794	-1.895 %
2.3	2163212	65.05169	-0.512 %
2.4	2551674	65.38679	-
2.5	2939390	65.4274	0.062 %

---

 | Show Table

DownLoad: CSV

Figure 9.  Optimum grid generation considered for each physics. "Pressure acoustics" physics: (a) whole geometry, (b) sliced magnified model. "Bio-heat transfer" physics: (c) whole geometry, (d) sliced magnified model.

DownLoad: Full-Size Img PowerPoint

3.   Validation

In order to verify the accuracy of our simulation methods, we proceed with the validation of the acoustic wave generation simulation, referred to as the "pressure acoustics" physics. Following that, we conduct the validation of the heat transfer simulation, known as the "bio-heat transfer" physics. These two validation steps, in addition to the comparison between PML and sound hard boundary conditions, are presented in the following sections.

3.1.   PML and sound hard boundary

The transducer was initially modeled using a sound hard boundary condition, which implies that the normal components of particle velocity and acceleration are zero at the boundary except at the transducer surface, the active element. However, recognizing the potential for unnecessary reflections with sound hard boundaries, we performed additional simulations using PML in Comsol Multi-physics. The PML effectively absorbed outgoing waves and minimized reflections. The comparison between PML and sound hard boundaries showed minimal differences in critical outcomes, such as maximum temperature and total acoustic pressure, confirming that our original results remain reliable. Details of this comparative analysis are provided in Table 3.

Table 3.  Comparison between PML and sound hard boundary conditions in FEM simulations results.

Parameters	Hard Boundary	PML
	(% error)
Maximum total acoustic	2.0026E7	1.9828E7
pressure (Pa)	(0.01 %)
Maximum sound intensity	1.1678E8	1.157E8
magnitude (Wm-2)	(0.0094 %)
Maximum Temperature (℃)	49.95	49.45
	(0.011 %)

---

 | Show Table

DownLoad: CSV

3.2.   Pressure acoustics phenomenon (Montienthong & Rattanadecho study ^[3])

To validate the initial stage of our simulation, we compare our results with those obtained by Montienthong & Rattanadecho ^[3], who conducted a 2D axisymmetric simulation of FUS in a cylindrical-shaped breast model, different from the ARBP model used in our study. Their model includes three tissue types—fat, glands, and muscle—within a larger water-filled cylinder, with a bowl-shaped transducer and two ultrasound frequencies (1 and 1.5 MHz) and tumor sizes (5 and 10 mm) considered. To ensure consistency, we adopted the same ICs and BCs, specifically simulating the 1 MHz frequency and 5 mm tumor scenario. We then compared the total acoustic pressure field and sound intensity magnitude distribution, focusing on the maximum values across the domain. The comparison between Montienthong & Rattanadecho's results and ours is presented in Figure 10.

Figure 10.  Comparison of the maximum total acoustic pressure and the maximum sound intensity magnitude obtained from our simulations (FEM and PINN) and the Montienthong & Rattanadecho's simulation.

DownLoad: Full-Size Img PowerPoint

It is apparent that there exist minor discrepancies between our simulation results and those of Montienthong & Rattanadecho. However, these small deviations can be considered acceptable within a tolerable margin of error.

3.3.   Heat transfer phenomenon (Montienthong & Rattanadecho study ^[3])

To verify the accuracy of our bio-heat transfer predictions, we compare our results with the simulations conducted by Montienthong & Rattanadecho ^[3], using the same geometric configuration and operating conditions. The temperature distributions predicted at the tumor's central region (focal point) during the ablation process are shown in Figure 11, alongside their results, demonstrating a good alignment between the two studies. Table 4 presents a comprehensive comparative assessment through relative error measurements, focusing on the highest overall acoustic pressure, the maximum sound intensity magnitude, and the maximum temperature achieved at the focal point during ablation. The relative error analysis reveals a satisfactory level of agreement between our simulations and those of Montienthong & Rattanadecho, as displayed in Table 4.

Figure 11.  Comparison of FEM and deep learning simulations results, and Montienthong & Rattanadecho simulation plots of the temperature distribution in the focal point during the time.

DownLoad: Full-Size Img PowerPoint

Table 4.  Comparison between our simulations results and Montienthong & Rattanadecho's simulation results ^[3].

			Montienthong &
Parameters	FEM	PINN	Rattanadecho's
	(% error)	(% error)	result
Maximum total acoustic	2.0026E7	1.9823E7	2.0228E7
pressure (Pa)	(-2 %)	(-2 %)
Maximum sound intensity	1.1678E8	1.156E8	1.1796E8
magnitude (Wm-2)	(-1 %)	(-2 %)
Maximum Temperature (℃)	49.95	50	49.9
	(0.1 %)	(0.2 %)

---

 | Show Table

DownLoad: CSV

4.   Results and discussion

In this section, we present the outcomes obtained from various simulated physical processes across different subsections. Initially, we showcase the field of total acoustic pressure and the distribution of sound intensity magnitude generated within the domain. Next, we illustrate the distribution of temperature in the breast model as a result of absorbing this acoustic energy. Finally, we present the extent of necrotized tissue and evaluate the amount of tumor and breast tissue that undergo necrosis due to the increase in temperature.

4.1.   Acoustic field generation

As previously mentioned, an attempt is made to generate a concentrated ultrasonic wave at the tumor site using a bowl-shaped acoustic transducer operating at a frequency of 1.1 MHz. Figure 12 displays the contours of total acoustic pressure and sound intensity magnitude, respectively. These contours depict the three cut planes (1 to 3) which include the x-y, z-y, and x-z surfaces passing through the center of the tumor. To provide a clearer depiction of the acoustic physics simulation outcomes, Figures 13 and 14 show the distribution of acoustic pressure and intensity magnitude along various x, y, and z directions, with the z-direction corresponding to the acoustic wave's propagation. As shown in Figure 13, the pressure amplitude is significantly higher at the tumor site due to focused acoustic waves, resulting in greater sound intensity and higher energy in the tumor region. Figure 14 highlights a prominent peak in sound intensity within the tumor, leading to rapid heat generation from dissipating this intense acoustic power. After simulating the pressure acoustics, the bio-heat transfer physics should be simulated to quantify the heat generated and the resulting temperature rise in the tumor region. The outcomes of this simulation are presented in the next subsection.

Figure 12.  Total acoustic pressure field, (a) FEM cut plane 1, (b) FEM cut plane 2, (c) FEM cut plane 3, (d) PINN cut plane 1, (e) PINN cut plane 2, (f) PINN cut plane 3. Sound intensity magnitude (Wm^-2), (g) FEM cut plane 1, (h) FEM cut plane 2, (i) FEM cut plane 3, (j) PINN cut plane 1, (k) PINN cut plane 2, (l) PINN cut plane 3.

DownLoad: Full-Size Img PowerPoint

Figure 13.  Absolute pressure distribution along lines passing through the tumor's center in different directions, (a) x-direction, (b) y-direction, (c) z-direction.

DownLoad: Full-Size Img PowerPoint

Figure 14.  Sound intensity magnitude distribution along lines passing through the tumor's center in different directions, (a) x-direction, (b) y-direction, (c) z-direction.

DownLoad: Full-Size Img PowerPoint

4.2.   Temperature distribution

The focused ultrasonic wave generates heat in the targeted area, raising the temperature, as shown in Figure 15(a)–(f) after 1 second on various cut planes. To better visualize the temperature rise in the tumor, Figure 16 presents temperature plots along lines passing through the tumor's center in different directions. By comparing Figure 15(a)–(f) with Figure 12(g)–(l), we observe that areas with higher sound intensity experience a greater temperature increase, with the tumor area reaching approximately 38.5 ℃ after a 1-second exposure. However, this temperature is too low for rapid tissue destruction, as thermal ablation, which causes tissue coagulation or necrosis, requires temperatures between 46 and 56 ºC sustained for a sufficient duration ^[54]. To visualize the necrotized area, the ablation time was extended to 10 seconds with a 10 nm transducer displacement, and the Arrhenius kinetic equation (Eq (2.6)) was used to calculate tissue necrosis during this period. The results of this analysis are presented in the subsequent subsection.

Figure 15.  Distribution of temperature after 1 second, (a) FEM cut plane 1, (b) FEM cut plane 2, (c) FEM cut plane 3, (d) PINN cut plane 1, (e) PINN cut plane 2, (f) PINN cut plane 3. Distribution of necrotic tissue fraction after 10 s of the ablation process, (g) FEM cut plane 1, (h) FEM cut plane 2, (i) FEM cut plane 3, (j) PINN cut plane 1, (k) PINN cut plane 2, (l) PINN cut plane 3.

DownLoad: Full-Size Img PowerPoint

Figure 16.  Temperature distribution after 1 s along lines passing through the tumor's center in different directions, (a) x-direction, (b) y-direction, (c) z-direction.

DownLoad: Full-Size Img PowerPoint

4.3.   Distribution of necrotized tissue fraction

Figure 15(g)–(l) illustrates the distribution of necrotic tissue after 10 seconds, showing the contour of the necrotic tissue fraction on different cut planes of the breast. As depicted, nearly the entire tumor is eradicated, with some surrounding normal tissue also affected. The maximum necrotic tissue fraction indicates complete necrosis in those areas. Comparing Figure 15(g)–(l) with Figure 15(a)–(f), a clear relationship emerges between the extent of necrosis and the temperature increase—regions with higher temperatures experience greater tissue ablation after 10 seconds of exposure. To provide further insight, Figure 17 presents plots showing the distribution of necrotic tissue fraction along lines passing through the tumor's center in various x, y, and z directions after 10 seconds of the ablation process.

Figure 17.  Fraction of necrotic tissue after 10 s ablation process along lines passing through tumor's center in different directions, (a) x-direction, (b) y-direction, (c) z-direction.

DownLoad: Full-Size Img PowerPoint

As shown in Figure 17, the proportion of necrotized tissue at the tumor's center reaches its maximum value of one but decreases as one moves away, indicating that not all tumor areas undergo complete ablation. To achieve full ablation, clinicians can adjust various parameters: briefly altering the acoustic pulse generator's position to change the focal point, increasing input pulse pressure to raise ablation temperature, or extending the duration of the acoustic pulse to expose the tumor longer. Multiple FUS ablation sessions can also target untreated tumor areas. Figure 18 presents plots of temperature and necrotic tissue fraction at the tumor's center throughout the 100-second ablation process, offering a comprehensive view of the procedure.

Figure 18.  Temperature and fraction of necrotic tissue at the tumor's center during ablation time: (a) Temperature change, (b) fraction of necrotic tissue.

DownLoad: Full-Size Img PowerPoint

As observed in Figure 18(a), during the initial 10 seconds of the existence of the acoustic pressure field, the temperature undergoes a rapid increase at the center of the tumor, elevating from 37 ℃ to approximately 70 ℃. Subsequently, upon cessation of the FUS source, the temperature at the tumor's center gradually decreases to 37 ℃, attributed to the heat transfer and cooling effect of the blood perfusion process. Conversely, depicted in Figure 18(b), the fraction of necrotic tissue at the tumor's center continually rises, reaching approximately 1 after 4 seconds. This indicates that approximately 100% of the tumorous tissues are ablated at the tumor's center.

4.4.   Realistic clinical approach

The previously described approach of focusing FUS solely at the tumor's center for 10 seconds followed by 90 seconds of thermal diffusion is inefficient for complete ablation, especially in larger tumors, and differs from current clinical practices. This method was designed to demonstrate the FUS process and our deep learning platform's capabilities. In this section, we adopt a more realistic clinical approach, showing how our platform can predict FUS ablation outcomes in actual clinical settings. In current clinical systems, a more effective method involves rapidly generating multiple lesions to cover the entire tumor volume by adjusting the transducer focus, thereby mitigating uncertainties related to heat diffusion and the heat sink effect caused by perfusion. In this study, we used a FUS therapy process with three faster lesions: one at the tumor's center and two others 3 mm away on the right and left sides along the x-axis. The transducer was moved and its focal point adjusted accordingly, with FUS waves generated for 6 seconds at each location, followed by 5-second intervals for repositioning. After treating the tumor from these three points, the simulation continued without any acoustic source until reaching 100 seconds.

The diagram in Figure 19(a)–(c) shows the distribution of sound intensity magnitude within the cut plane 1 for FEM and (d to f) for PINN, with the x-z surface passing through the tumor's midpoint for each of the three lesion processes occurring between 0–6, 11–17, and 22–28 seconds. Additionally, Figure 19(g)–(r) presents the distribution of temperature and the fraction of necrotic tissue at the conclusion of each lesion process within cut plane 1 for FEM and PINN, concluding at 6, 17, and 28 seconds for the first, second, and third lesions, respectively. These processes involve three distinct focal points where a significant sound intensity magnitude of approximately 4.5E+6 Wm^-2 is generated, leading to temperature increases to about 63, 65, and 64 ºC at different points within the tumor, contributing to a more uniform necrosis of the tumor tissue.

Figure 19.  (a)–(f) sound intensity magnitude in the cut plane 1 during each of the three lesion processes for FEM and PINN, (g)–(l) temperature distribution in the cut plane at the end of each lesion process for FEM and PINN, (m)–(r) fraction of necrotic tissue distribution in the cut plane 1 at the end of each lesion process for FEM and PINN.

DownLoad: Full-Size Img PowerPoint

The figure demonstrates that after 28 seconds of FUS therapy, the fraction of necrotic tissue reaches 1 in nearly all regions of the tumor, indicating complete ablation. To better understand these results and evaluate the impact of this clinically realistic approach, temperature and necrotic tissue fraction are plotted along a line passing through the tumor's center in the x direction (Figure 20) at various time points during the ablation process, including the end of each lesion process and the entire therapy. Figure 20(a) shows the temperature increase at different points within the tumor—center, right, and left—resulting in a more evenly distributed high-temperature pattern. Comparing Figure 20(b) with Figure 17(a), it is evident that this approach achieves a fraction of necrotic tissue reaching 1 across almost all regions, unlike the previous method, where this fraction did not reach 1 at the tumor's edges, highlighting the new approach's effectiveness in uniformly and completely eradicating tumor tissues.

Figure 20.  Temperature and the fraction of necrotic tissue after 6, 17, 28 and 100 s along the line passing through the tumor's center in the x direction.

DownLoad: Full-Size Img PowerPoint

Increasing the number of lesion points across different parts of the tumor leads to more favorable outcomes without requiring extended lesion durations, aligning with current clinical practices. This deep learning platform can simulate various FUS ablation processes, accommodating different lesion numbers, durations, locations, and acoustic settings, making it a valuable tool for predicting diverse FUS therapy results and selecting the most efficient procedure based on tumor size and location. Notably, 3D simulations offer significantly more precise outcomes than 2D simulations due to the increased complexity and incorporation of three spatial dimensions, resulting in higher sound attenuation and a more pronounced temperature increase within the tissue.

The augmented temperature increase in our simulation, evident when comparing our results with those of Montienthong & Rattanadecho's study ^[3], highlights the advantages of 3D simulations. While our study shows a temperature rise to approximately 64–65 ℃ (Figure 20(a)), their research only reached 50–51 ℃ despite using the same acoustic frequency and similar focal lengths (Figures 14 and 15 in their paper). The 3D simulations allow for incorporating the intricate and realistic geometry of the breast model, including asymmetric details that significantly affect acoustic wave propagation and heat transfer, leading to more accurate results than 2D axisymmetric models. The PINN model's ability to generalize to different tumor types and breast anatomies is crucial for clinical applicability. Although our study demonstrates the model's effectiveness for a specific anatomy and tumor type, generalization to other scenarios requires defining new anatomical geometries and FUS specifications. The training dataset was based on specific configurations of an anatomical phantom, capturing particular scenarios of acoustic pressure, sound intensity, and temperature distributions. However, this dataset doesn't cover the full range of breast anatomies and tumor types, meaning the model may not perform optimally for cases significantly different from those represented in the training data, requiring precise input configurations for accurate FUS treatment predictions.

When scaling the PINN approach to more complex anatomical structures, several challenges arise, including increased model complexity and computational demands. However, the modular nature of PINNs makes them well-suited for handling such complexity. One approach to managing complexity is hierarchical modeling, where the anatomical structure is divided into smaller, more manageable subdomains. Each subdomain can be modeled individually using the PINN framework, and the results can be integrated to form a comprehensive model of the entire structure. Adaptive mesh refinement techniques can be used to focus computational resources on regions with high complexity or significant changes, such as tumor boundaries or areas with steep temperature gradients. This allows for more detailed modeling in critical areas without overwhelming the computational resources.

4.5.   Validation and computational efficiency

To validate the accuracy of the proposed PINN model, we compared its predictions against those from FEM simulations. The validation process was carried out across different clinical scenarios, varying key parameters such as tissue properties, FUS power, and blood perfusion rates. We assessed the model's performance using metrics such as root mean squared error (RMSE), mean absolute error (MAE), and the correlation coefficient (R²). The following table (Table 5) summarizes the validation metrics in various scenarios.

Table 5.  Validation metrics for comparison with FEM results in various scenarios.

Scenario	RMSE (℃)	MAE (℃)	Correlation coeddicient (R2)
Baseline (nominal tissue properties, FUS power)	0.05	0.03	0.98
±10% Variation in Thermal Conductivity	0.07	0.05	0.97
±10% Variation in Specific Heat Capacity	0.06	0.04	0.97
±10% Variation in FUS Power	0.08	0.05	0.96
±10% Variation in Blood Perfusion Rate	0.06	0.04	0.98

---

 | Show Table

DownLoad: CSV

One of the key advantages of the PINN model is its computational efficiency compared to conventional FEM techniques. A detailed comparison of computational time and resource usage between the PINN model and FEM simulations, as shown in Table 6, demonstrates that the PINN model significantly reduces both computational time and memory usage. Specifically, the PINN model completed the simulations in just 5 hours, whereas FEM simulations required approximately 115 hours. Additionally, the PINN model used 8 GB of memory, compared to 16 GB for FEM simulations. These substantial differences underscore the efficiency of the PINN model, making it particularly advantageous for clinical applications where timely results are essential. By leveraging this computational efficiency, the PINN model enables rapid and accurate simulations of FUS ablation, making it a valuable tool for treatment planning and optimization.

Table 6.  Computational time and resource usage comparison.

Method	Computational time (hours)	Memory usage (GB)
PINN Model	5	8
FEM Simulation	115	16

---

 | Show Table

DownLoad: CSV

4.6.   Uncertainty quantification

Uncertainty quantification (UQ) is crucial for validating and ensuring the robustness of predictive models, especially in clinical applications. In this study, we addressed uncertainties in model parameters and input data to enhance the reliability of PINN predictions. Parameters such as tissue properties (e.g., thermal conductivity, specific heat capacity) and blood perfusion rates, as well as variations in input data like FUS wave intensity and focus, can significantly impact model predictions. To account for these uncertainties, we conducted Monte Carlo simulations, running multiple simulations with parameters sampled from their respective probability distributions. The results provided statistical measures such as mean, variance, and confidence intervals for predicted quantities like temperature distribution. Additionally, sensitivity analysis was performed to identify the parameters most significantly affecting model outputs, with findings summarized in Table 7. This involved varying one parameter at a time while holding others constant to assess how each affects predicted temperature and necrotic tissue distribution. Bayesian methods were also used to update parameter probability distributions based on prior knowledge and observed data, incorporating prior uncertainties and improving the model's predictive accuracy.

Table 7.  Sensitivity analysis results.

Parameter	Range	% Change in maximum temperature	% Change in necrotic tissue fraction
Thermal conductivity	0.4–0.6 W/m·K	±5%	±3%
Specific heat capacity	3500–4000 J/kg·K	±4%	±2%
Blood perfusion rate	0.001–0.005 s-1	±7%	±5%
Initial tumor temperature	36.5–37.5℃	±2%	±1%
FUS power	30–50 W	±10%	±8%

---

 | Show Table

DownLoad: CSV

The uncertainty quantification results are presented in Table 8, which includes the mean and standard deviation of key predicted quantities, as well as the 95% confidence intervals. These results demonstrate the robustness of the PINN model predictions despite the presence of uncertainties in the model parameters and input data. The narrow confidence intervals indicate that the model is relatively insensitive to parameter variations, which enhances confidence in its clinical applicability. For example, we explored the effect of a ±5% variation in thermal conductivity and specific heat capacity on the temperature distribution within the tumor. The simulations show that a 5% decrease in thermal conductivity leads to a 3% increase in peak tumor temperature, while a 5% increase in thermal conductivity results in a 2% reduction in peak temperature. This demonstrates that the tumor temperature is sensitive to variations in tissue properties, but the overall treatment outcome remains within clinically acceptable bounds, with full ablation achieved in all cases.

Table 8.  Uncertainty quantification results.

Quantity	Mean	Standard deviation	95% Confidence interval
Maximum temperature (℃)	70.0	2.5	[65.2, 74.8]
Necrotic tissue fraction (%)	90.0	5.0	[80.2, 99.8]

---

 | Show Table

DownLoad: CSV

Uncertainty in the applied FUS power was also analyzed, with a ±10% variation around the nominal value. A 10% increase in FUS power led to a 6% increase in the extent of necrotic tissue, while a 10% reduction in power resulted in a 4% decrease in necrotic tissue. While these variations could impact treatment efficiency, the model demonstrated that tumor ablation remained effective even at lower power levels, showcasing its robustness to fluctuations in FUS power during clinical procedures. We considered a range of blood perfusion rates from 0.001 to 0.005 s^-1 to simulate how variability in perfusion impacts the cooling effect during ablation. Higher perfusion rates led to a slight reduction in the maximum temperature (approximately 5%) and a delayed necrosis onset. However, the tumor region still reached the required thermal dose for complete ablation in all cases, further supporting the reliability of the model across different physiological conditions.

These analyses demonstrate that while uncertainty in model parameters can lead to variations in the exact temperature distribution and necrotic tissue extent, the overall therapeutic outcomes—namely, successful tumor ablation—are robust. This highlights the model's ability to handle clinical uncertainties and maintain treatment effectiveness. Moving forward, this uncertainty quantification will aid in establishing safety margins for clinical treatments and provide confidence in the model's predictions when used in real-world FUS therapy.

4.7.   Limitations and future work

While our study provides a comprehensive simulation of FUS ablation using a PINN model, it is important to acknowledge certain limitations, particularly the use of the traditional Pennes bio-heat equation, which assumes homogeneous or isotropic blood perfusion and may not accurately reflect the heterogeneous and anisotropic conditions within breast tissue. Recent studies have suggested modifications to account for spatial variability in blood perfusion ^[55,56], especially in tumor hypoxia, where low blood perfusion can affect thermal energy deposition and reduce sensitivity to ablation treatments. Moreover, the current model assumes linear acoustic propagation, which may not accurately capture the complexities of FUS, where nonlinear effects can significantly influence the pressure distribution and temperature rise at the focus, potentially leading to discrepancies between the simulated and actual acoustic pressure values, particularly in high-intensity regions.

Additionally, it is important to consider the continuous regeneration of living tissues due to oxygen supply through arterial blood, which counterbalances thermal degradation under quasi-static thermal conditions. Heating of cellular and biological tissues can lead to reversible (repairable) or irreversible (lethal) thermal cell death. At the interface between tumor and healthy tissue, the regeneration of healthy cells triggers an immune response that suppresses, prevents, and restricts further thermal damage within the damage bounds of $\Omega \leq$ 1. This partial self-regeneration of normal tissues at the tumor periphery, facilitated by continuous oxygen supply, should not be overlooked when modeling thermal damage kinetics.

The current bio-heat model does not consider the effects of large blood vessels and thermal relaxation in bio-tissues, both of which are significant in heat transfer. Large blood vessels can act as heat sinks, rapidly dissipating heat and affecting thermal distribution, while thermal relaxation accounts for the time delay in heat conduction. Ignoring these factors can lead to inaccuracies. Future work should integrate a modified bio-heat equation that considers heterogeneous blood perfusion, large blood vessels, and thermal relaxation into our PINN framework. This enhancement would improve the accuracy and clinical relevance of our simulations, allowing for a more detailed analysis of thermal energy distribution and its effects on tumor ablation, providing valuable insights for optimizing FUS treatment protocols.

While the current model demonstrates effective generalization within the scope of the training data, future work will focus on expanding the dataset to improve its applicability across a broader range of breast anatomies and tumor types. This will be achieved through the inclusion of additional 3D models from datasets such as UWCEM, along with patient-specific MRI-derived models representing different breast densities, tumor sizes, and locations. Furthermore, we plan to apply transfer learning techniques to enhance the model's adaptability to new clinical cases. Data augmentation strategies, such as simulating variations in tissue properties and tumor characteristics, will also be employed to further enrich the dataset. Lastly, rigorous testing on diverse clinical cases will validate the model's generalization capability and robustness in real-world scenarios.

Our PINN-based solution operates efficiently in a mesh-free domain, providing a simplified framework that offers preliminary insights into the FUS process. However, we recognize that this model does not fully capture the complexities of real-world FUS therapy scenarios. Future work will focus on enhancing the model to incorporate more detailed tissue properties, patient-specific anatomy, and dynamic physiological conditions to improve its applicability to clinical settings. By addressing these limitations and incorporating more advanced modeling techniques, we aim to improve the predictive power and applicability of our deep learning platform for FUS ablation therapies and extend the model to other types of cancer.

5.   Conclusions

This study presents a comprehensive deep learning framework based on PINN for 3D simulation of FUS ablation in breast tissue, accurately replicating key aspects such as acoustic pressure, sound intensity, temperature distribution, and necrotic tissue. We validated our model against existing numerical findings and FEM simulations, demonstrating its effectiveness in simulating FUS ablation with an anatomically accurate breast phantom. The results show that a bowl-shaped acoustic transducer can effectively focus ultrasound waves to raise the temperature at the tumor site, leading to significant tumor necrosis.

We explored a clinically relevant approach involving multiple lesions, achieving a more uniform temperature distribution and effective tumor ablation, demonstrating the advantages of 3D simulations over 2D models in capturing FUS ablation processes. Our deep learning model enables personalized treatment planning, potentially improving clinical outcomes by providing precise predictions of temperature elevation and tumor ablation. However, the model's effectiveness depends on the quality of training data and its generalizability to diverse patient cases, requiring extensive validation and clinical trials before implementation.

Future research should focus on expanding the dataset, improving model interpretability, and integrating the model into clinical workflows. Continuous monitoring and adaptation of the model will ensure its long-term utility in clinical settings.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by the Korea-Canada Artificial Intelligence Joint Research Center at the Korea Electrotechnology Research Institute (Operation Project: No.295 22A03009/22A03007), which are supported by Changwon and Busan City, Korea.

Conflict of interest

The authors declare there is no conflict of interest.