Mesh-Free Modeling of Heat Transfer Dynamics for Rapid Assessment of Necrosis Zones in Hepatic Tumor Radiofrequency Ablation Using Physics-Informed Neural Networks

Hepatocellular carcinoma and secondary liver metastases are particularly unsuitable for surgical resection due to the location of the tumor or underlying cirrhosis [1]. Radiofrequency ablation (RFA) has emerged as a primary and minimally invasive thermal therapy that has become the standard of care for unresectable liver tumors.

RFA induces coagulative necrosis by converting electromagnetic energy into heat through the percutaneous insertion of a needle electrode directly into the target tissue. Alternating current operating between 460 and 500 kHz induces rapid ionic agitation and molecular friction in adjacent tissue [2] to generate continuous heat to raise the local tissue temperature above 50 °C for a period that ensures irreversible protein denaturation and coagulative necrosis of malignant tumor cells [3].

The efficacy of RFA in liver tissue is fundamentally governed by the dynamics of complex spatiotemporal heat transfer. The mathematical description of heat transfer during RFA is mainly based on the Pennes bioheat equation [4]. This partial differential equation (PDE) combines standard thermal conductivity with the convective cooling effect of blood flow. The liver is characterized by a dense and highly active microvascular network. During thermal treatments, the mechanism by which local blood perfusion provides a strong convective cooling is widely documented in the literature as a heat sink effect [5]. As the electrode generates heat, nearby blood vessels actively dissipate thermal energy from the targeted ablation site. This dynamic cooling prevents the tissue from reaching necrosis temperature at the boundaries of the tumor, resulting in incomplete ablation and subsequent recurrence of the local tumor. Recent clinical literature has emphasized that the recurrence of post-treatment disease remains a significant limitation [6], [7]. Comparative ex vivo models emphasize that monopolar RFA is particularly vulnerable to this heat-sink phenomenon compared to other modalities, significantly distorting the targeted ablation volume and density [8].

Analytical solution is almost impossible due to the nonlinear nature of tissue parameters for realistic 2D or 3D hepatic geometries and the complex spatial distribution of radiofrequency energy source [9]. Recent analyses have physically substantiated this recurrence of risk by revealing that large adjacent vessels not only restrict lesion expansion but also shift the thermal hotspot and cause dangerous underdose margins [10]. Achieving complete tumor treatment transforms from a static geometric planning issue into a dynamic control problem, further complicated by sudden changes of tissue property like impedance spikes due to localized vaporization [5], [11], [12]. Consequently, accurate pretreatment modeling of these spatial and temporal thermal dynamics is critical for treatment success [13]. Within this clinical context, the integration of computer modeling into treatment optimization has emerged as a rapidly growing research trajectory, shifting the paradigm from purely experiential intuition to measurable and patient-specific targets [9], [14].As illustrated in Figure 1, the contraction of the ablation zone due to the heat sink effect is a significant factor that affects treatment outcomes [15].

The estimation of temperature distribution in biological tissues has relied on network-based numerical techniques such as finite element method (FEM) and finite difference method (FDM) [16]. Computational simulations function as essential components of advanced clinical Decision-Making Systems [17]. Advanced simulations enable the optimization of critical parameters such as electrode placement, applied power, and exposure time during the planning phase, thereby improving therapeutic outcomes and minimizing damage to surrounding healthy tissues [16].

The FEM and FDM methods offer high physical fidelity, but they require extensive spatial decomposition and extremely small timesteps to meet the stability criteria. This constraint requires long computational times which often exceed the duration of the clinical procedure, and thus the computational method becomes impractical for the intraoperative decision support. A recent systematic review of liver ablation simulations confirms this bottleneck, noting that while classical FEM models are structurally preferred for solving coupled bioheat equations in complex geometries, they suffer from high computational costs and a persistent standardization gap across validation metrics [18]. To address delays in these clinical workflows, recent parallel efforts have introduced graphics processing unit (GPU)-accelerated pre-computed frameworks [16] and adaptive algorithms using temperature-impedance feedback [19].

To shorten the computation times, recent studies have applied fully data-driven machine learning algorithms to predict temperature distributions in biological tissues. Ren et al. [20] trained a two-branch artificial neural network model with data obtained from a computational model and tested the model by data from ex vivo experiments. They reported that the model efficiently predicted lesion size RFA. The two-branch approach was argued to have improved the ability of their model to fit complex relationships through activation functions and nonlinear combination features [20]. However, standard neural networks operate as black-box models. As they rely entirely on large and pre-labeled datasets that infer relationships between variables and intrinsically disregard the fundamental laws of thermodynamics [21], they are prone to producing physically ill-posed solutions and biologically impossible thermal gradients [1], [22]. In biomedical, where experimental data is often sparse, noisy, or difficult to obtain in vivo, these purely data-driven models have drawbacks in generalizing outside their training areas. Generating physically inconsistent estimates poses serious risks to clinical practice, such as overestimating or underestimating the ablation site directly impacting a patient’s safety.

To address the limitations of both classical numerical methods, solvers and standard data-driven machine learning, the 2D Physics-Informed Neural Networks (PINNs) framework was introduced to solve the spatiotemporal Pennes bioheat equation for hepatic RFA simulation. As a robust hybrid model that combined the high physical fidelity with the speed of neural networks, PINNs bypassed traditional stability constraints by continuously evaluating spatial and temporal coordinates (x, y, and t) without relying on a discrete grid. They embedded the governing PDEs, along with the corresponding initial and boundary conditions, directly into the neural network’s Loss Function [23] and forced the algorithm to minimize the residual of physical laws during the training phase to ensure that the outputs adhered to thermodynamic principles [24].

As summarized in Table 1, the PINN framework eliminates the mesh-dependency and time-stepping constraints of classical methods and facilitates real-time prediction of the ablation zone. Recent advancements in digital twin frameworks and patient-specific validation have emphasized that biomarkers such as tissue composition, steatosis, and complex vascular structures significantly alter ablation volumes, pushing the field toward highly personalized “margin and safety” assessments [9], [25], [26]. However, incorporating these highly heterogeneous and patient-specific boundaries into traditional FEM requires extremely labor-intensive and computationally heavy mesh generation and refinement. By eliminating mesh-dependency entirely, the proposed PINN framework seamlessly accommodates spatially varying biophysical parameters, such as the localized Gaussian perfusion profiles of large vessels, without the geometric constraints or remeshing bottlenecks of classical methods. Ultimately, the integration of PINNs into clinical workflows could significantly improve the accuracy of treatment planning and allow clinicians to optimize electrode placement and power settings based on patient-specific perfusion characteristics [24].

The governing equation for the balance of thermal energy in biological tissues is the Pennes bioheat equation that combines conductive heat transfer with the convective cooling effects of blood perfusion. For a two-dimensional spatial domain with time dependence, the formulation is expressed in Eq. (1):

where, $\rho$ is the density of the tissue, $c$ is the specific heat capacity, and $k$ is the thermal conductivity. The term $\omega_b \rho_b c_b\left(T_a-T\right)$ defines the cooling effect caused by blood perfusion, which is significant in highly vascularized organs such as the liver. In this term, $\omega_b$ is the blood perfusion rate, $\rho_b$ is the density of the blood, $c_b$ is the specific heat of the blood, and $T_a$ is the arterial blood temperature. The variables $Q_m$ and $Q_{R F A}$ represent the metabolic heat generation and the external heat source from the RFA electrode, respectively.

Accurate simulation requires appropriate thermophysical and physiological parameters for human hepatic tissue. The liver exhibits strong heat sink characteristics due to extensive microvascular networks; therefore, the perfusion term dominates the thermal response during ablation. Standard literature values, detailed in Table 2, were used for the baseline model.

The heat generated by the radiofrequency electrode is modeled as a volumetric spatial heat source. In practical RFA applications, the dissipation of electromagnetic power decreases rapidly with distance from the tip of the electrode, as formulated in a continuous and differentiable form below:

where, $P_{\max }$ is the maximum applied power density at the electrode center, and $\sigma^2$ is a variance parameter controlling the spatial spread of the thermal energy.

The computational domain $\Omega$ is defined as a two-dimensional square cross-section of hepatic tissue, measuring [-0.05, 0.05] meters in both the $x$ and $y$ directions, the geometric center, or which correspond to the placement of the RFA electrode.

The initial condition assumes uniform thermal equilibrium at arterial body temperature prior to the onset of ablation at $t=0$ :

Dirichlet boundary conditions are applied at the boundaries of the domain $(\partial \Omega)$, under the assumption that the boundaries are sufficiently far from the heat source to remain unaffected by the ablation process:

An explicit Forward-Time Central-Space (FTCS) finite difference scheme was implemented to compute the distribution of baseline temperatures. The computational domain is discretized into a uniform twodimensional grid with spatial steps $\Delta x$ and $\Delta y$ and a temporal step $\Delta t$. Spatial derivatives, on the other hand, are approximated using second-order central differences, while the time derivative is a first-order forward difference. The discrete update equation for $T_{i, j}^n$, the temperature at spatial node $i, j$ at time step $n$ is defined as in Eq. (5):

where, $\alpha=k /(\rho c)$ is the thermal diffusivity of the hepatic tissue.

To ensure numerical stability in the explicit FTCS scheme, the time step $\Delta t$ is strictly constrained by the Courant-Friedrichs-Lewy criterion for the two-dimensional conduction problem:

This requirement of stability restricts $\Delta t$ to small values, which consequently increases the computational time required to simulate a standard 10 minute ablation procedure.

The predictive accuracy and convergence stability of Physics-Informed Neural Networks (PINN) are highly sensitive to the capacity of the network architecture, i.e., the number of hidden layers and neurons, and the learning rate in the optimization process. In order to model highly non-linear biophysical constraints such as interaction of latent heat, impedance roll-off mechanism, and vascular heat sink effect in the RFA process, it is critical for the network to have sufficient representation capacity.

However, an excessively large network architecture brings about the risks of overfitting and high computational costs, while an inadequate architecture cannot solve complex thermal gradients due to underfitting [8] ( Figure 2). Therefore, a systematic hyperparameter sensitivity analysis was carried out on three different hidden-layer topologies and three different initial learning rates to determine the optimum network configuration. These nine different scenarios were pre-trained with the standard Adam optimization with 3 000 steps and the loss values and computation times were recorded.

The heatmap of final loss and detailed results are presented in Figure 3 and Table 3, respectively.

The training efficiency and convergence stability of different network depths/widths and learning rates in Table 3 show that all tested architectures converged satisfactorily and that the established mathematical model and the physical constraints were stable across all network architectures.

Obviously, the ability to minimize the PDE residue increased consistently with the number of layers and neurons. The lowest training loss was obtained in the 5 $\times$ 60 configuration with a learning rate of $10^{-3}$ whereas the 4 $\times$ 50 network had reached almost the same final loss within nearly half of the time. Therefore, the 4 $\times$ 50 configuration was adopted for the simulations and a fully-connected feedforward neural network, shown in Figure 4, was constructed to map the spatiotemporal coordinates to the tissue temperature, $T(x, y, t)$.

The network architecture consists of an input layer with three neurons that represent the independent variables ($x$, $y$, and $t$), multiple hidden layers to capture the nonlinear dynamics of the perfusion term, and a single neuron output layer for the predicted temperature, The hyperbolic tangent activation function, which is continuously differentiable, is used across all hidden layers to compute the spatial and temporal gradients of $T$.

The training of PINN is driven by the minimization of a composite loss function, $\mathcal{L}$ that penalizes deviations from the governing physical laws, the initial state, and the boundaries:

where, $\omega_{P D E}, \omega_{I C}$, and $\omega_{B C}$ are weighting coefficients used to balance the gradient contributions of each term.

The physical loss, $\mathcal{L}_{\text {PDE }}$, enforces the Pennes bioheat equation across a set of collocation points, $N_f$, sampled within the spatiotemporal domain, and is mathematically expressed as:

where, $\hat{T}$ is the temperature predicted by the neural network. The terms governing spatial gradients $\left(\nabla^2 \hat{T}\right)$ and temporal gradients $\left(\frac{\partial \hat{T}}{\partial t}\right)$ are computed exactly using automatic differentiation.

The initial condition loss, $\mathcal{L}_{I C}$, and the boundary condition loss, $\mathcal{L}_{B C}$, are calculated using mean squared error (MSE) over sets of points $N_{I C}$ and $N_{B C}$, sampled at $t$ = 0 and the domain boundaries, respectively:

By minimizing this composite loss, the network converges to a state that satisfies the PDE without requiring a pre-computed labeled dataset.

Standard analytical solutions to the Pennes equation often assume constant tissue properties and linear heat generation. However, clinical RFA involves complex non-linear phenomena that significantly affect thermal distribution. To enhance the physiological accuracy of the PINN framework, two critical biophysical constraints for temperature-dependent impedance roll-off and spatially heterogeneous perfusion were integrated directly into the governing equations.

The heating efficiency in biological tissues is physically limited by the phase change of intracellular water. As the tissue temperature approaches 100 °C, vaporization leads to desiccation, which causes a sharp increase in electrical impedance and effectively stops the current flow (the "roll-off" effect). To mimic this self-limiting behavior, a smooth sigmoid-based switching function, $\sigma(T)$, was integrated into heat source term, $Q_{R F A}$ without introducing nondifferentiable discontinuities that would destabilize the training phase of PINN:

The switching function is defined as:

where, $T_{\text {crit }}$ is the vaporization threshold of 100 °C and $k$ is the steepness control of transition. This mechanism automatically approaches zero as the local temperature exceeds $T_{\text {crit }}$ and prevents overheating.

The cooling effect of large blood vessels is the primary cause of incomplete ablation and recurrence of local tumor. To simulate the convective heat sink effect, the baseline blood perfusion rate $w_b$ was changed from a scalar constant to a spatially varying function $w_b(x, y)$. The presence of a large liver vessel was modeled using a Gaussian perfusion profile as follows:

where ($x_v, y_v$) are the vessel coordinates, $\delta$ is the radius of the vessel, and $\alpha$ is the intensity factor for the high flow rate within the vessel. By embedding this function in the PDE residual, PINN learns to adjust the thermal gradients specifically in the perivascular region and creates a more realistic deformation of the ablation zone.

Although the integrated impedance roll-off mechanism suppresses temperatures over 100 1111, the latent heat of vaporization has been incorporated into the PDE formulation.

Liver tissue does not evaporate at an exact point but spreads over a narrow temperature range as it contains dissolved salts and proteins and has differences in intracellular pressure.

In order to accurately model the enthalpy of vaporization without tracking a moving boundary, an effective heat capacity method that spreads the total latent heat over a sinusoidal range, as in eq. (14), should be used.

Since this function should be valid only in the phase-change range, using a piecewise function in PINN would disrupt the derivative continuity and collapse the training, a Gaussian Distribution was employed instead.

By doing so, the sinusoidal distribution was mathematically imitated almost exactly but an infinitely differentiable function was achieved:

where, $\sigma$ is the standard deviation parameter of the Gaussian distribution and determines the width of the temperature band in which this evaporation occurs. In the present PINN model, $\sigma=1$ was used, which allowed the latent heat to be absorbed in the range $\pm$ 1 °C.

The amplitude of the heat capacity peak was mathematically calibrated such that the integral of the excess capacity over the phase change range equaled the theoretical enthalpy of vaporization:

This formulation provides a thermodynamically consistent prediction of the ablation zone by ensuring that the exact amount of electromagnetic energy required to vaporize tissue water is consumed by the model, before any temperature rise above 100 °C.

The clinical endpoint of RFA is the induction of irreversible coagulative necrosis. This time-dependent biological damage is governed by the Arrhenius formulation [27]:

where, $\Omega$ is the degree of tissue damage, $A$ is the frequency factor, $E_a$ is the activation energy, $R$ is the universal gas constant, and $T$ is the temperature in Kelvin. Cell death occurs only when $\Omega \geq$ 1 and the percentage of necrosis increases with $\Omega$. A 99% cell death is guaranteed with values $\Omega \geq$ 4.6. In hepatic tissue, the threshold is typically crossed when temperatures exceed 50 °C for approximately one minute, or instantaneously at temperatures above 60 °C due to rapid protein denaturation [3], [28].

The computational framework for the PINN model was implemented in Python using the DeepXDE library. The spatial domain and time interval were defined, and collocation points were generated using Latin Hypercube Sampling to ensure uniform coverage across the 2D + t space ( Figure 5).

To rigorously capture the steep thermal gradients near the heat source singularity and prevent the nonphysical dampening of peak temperatures at the tip of the electrode, the uniform Latin Hypercube Sampling was supplemented with a spatially adaptive collocation point refinement strategy. A dense cluster of auxiliary training points was dynamically concentrated within a 5 mm radius of the electrode center ($x$ = 0, $y$ = 0). This localized over-sampling ensures that the neural network allocates sufficient representational capacity to resolve the singularity, hence strictly preserving the geometric fidelity of the expanding 50 °C isotherm which defines the ablation margin.

The training set then consisted of 8000 interior domain points, 2000 boundary points, and 2000 adaptive anchor points near the tip. The initial condition points were removed from the training set where they were treated as a penalty term in the composite loss function in order to avoid competing gradients and temporal error propagation. The absolute physical fidelity at $t$ = 0 was ensured by exact enforcement formulation.

By applying a spatiotemporal output transform ansatz:

The initial condition was satisfied by the network architecture itself. Since the integrated circuit (IC) error is mathematically fixed to zero, it did not require optimization.

A two-stage training strategy was used in the optimization process to handle the gradient differences between the spatial and temporal domains. Initially, the Adam optimizer was employed for 15000 iterations with a learning rate of $1 \times$ $10^{-3}$ to rapidly reduce the global loss and avoid local minima. This stage rapidly reduced the global loss to $4.50 \times$ $10^{-1}$ but struggled to achieve deeper convergence due to the stiffness of the bioheat equation.

The training process leveraged 8000 interior domain points and 4000 boundary and initial condition points.

The subsequent application of the second-order Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer for an additional 11 000 steps effectively finetuned the network weights and accelerated convergence near the exact solution. The training concluded in 3 348 seconds and achieved a final composite loss of $1.5 \times$ $10^{-1}$. The loss history curve presented in Figure 6 demonstrates a stable logarithmic decay in the governing PDE residual, which indicates that the neural network successfully embedded the thermodynamic constraints rather than merely memorizing boundary conditions.

The progression of the thermal field during the simulated 10 min ablation procedure is illustrated through spatiotemporal heatmaps at 60,300 , and 600 seconds. The 50 °C isothermal contour, representing the clinical boundary for irreversible coagulative necrosis, expands radially outward from the electrode center.

The spatial temperature profile along the $y$ = 0 section shows the Gaussian nature of thermal propagation. For the non-vascular heat sink scenario (Figure 7a), the central temperature reached above the necrosis threshold 50 °C even at $t$ = 60 seconds and approached a maximum of 54 °C at $t$ = 600 seconds. The thermal gradient shows a significant decrease when moving outward active electrode, primarily with the low thermal conductivity of the tissue and the convective heat absorption by liver blood perfusion.

To assess the capability of the model in the complex clinical field, a vascular heat sink effect was modeled by adding a spatially varying perfusion function that introduced a large vessel at $x$ = 1.5 cm near the electrode center. As can be seen in Figure 8b, the vessel induced a great variation of the thermal distribution. While the ablation zone expanded symmetrically in the $x$ $<$ 0 region, the proximity of the heat sink in the $x$ $>$ 0 region severely skewed the temperature field. High convective cooling induced by the vessel prevented the vicinity from reaching the cell necrosis temperature of 50 °C. At $t$ = 600 s, the ablation edge showed a concave indentation at a depth of about 5 mm. The contour plots show that the PINN framework correctly resolves directionally dependent heat sink effect without requiring explicit meshing or mesh refinement around the vessel.

To quantitatively validate the clinical safety of the PINN predictions and ensure that the functional dampening at the heat source singularity does not compromise the ablation margin, spatial overlap and boundary distance metrics were computed specifically for the critical 50 °C coagulative necrosis contour against an explicit 2D FDM baseline seen in Figure 9. The spatial agreement was evaluated using the Dice Similarity Coefficient (DSC) below:

The calculated DSC values of 0.9799 and 0.9041 showed that the PINN framework achieved a near-perfect volumetric overlap of the predicted necrosis zones. The maximum boundary deviation was measured via the Hausdorff Distance (HD) to be defined as:

The Hausdorff Distance for the 50 °C contour was restricted to 1 mm, and 3.162 mm, respectively for the cases modeled. As this minor boundary deviation falls below the standard clinical safety margin requirement, it shows that functional dampening at the extreme peak of the singularity at the tip of the electrode does not compromise or artificially shrink the predicted ablation boundary, hence confirming that the mesh-free continuous approximation by PINN could reliably resolve the therapeutic boundary.

The difference between the spherical ablation region in Figure 10a and the deformed shape observed in Figure 10c shows the critical importance of patient-specific modeling. The indentation in highly deformed necrosis temperature zone that is likely to have viable tumor residue would not be foreseen in traditional RFA planning that assumes homogeneous heat distribution; this would result in inadequate treatment and potential local recurrence. The ability of PINN to model the heat seat in real time offers significant advantages for preoperative planning and intraoperative decision making in adjusting electrode placement or power protocols.

Analysis of transient temperature evolution presented in Figure 11 shows the variation of temperature at the center point of the electrode ($x$ = 0, $y$ = 0) and at the periphery of 1 cm in $x$- and $y$-axes. The central point experienced a fast and nonlinear heating while reaching the necrosis threshold soon after 100 seconds. The peripheral tissue at a 1 cm distance exhibited a delayed heating response and crossed the necrosis threshold later in the procedure, due to the spatial decay of the electromagnetic energy and continuous cooling from blood perfusion.

In the first scenario where no vascular heat sink effect was induced, the peripheral tissue temperature remained over the necrosis threshold for about twice the necessary duration of 1 minute. In the presence of a large vessel near the tumor, however, the decline in the temperature of the peripheral tissue was dramatic and remained totally under the necrosis threshold.

It is important to note that predicted peak temperatures in models represent theoretical values derived from the bioheat equation. In clinical practice, tissue temperature was physically the phase change of intracellular water. When it reached this boiling point, tissue vaporization caused a sharp increase in electrical impedance (the “roll-off” phenomenon), which acted as an electrical insulator and limited further heating [29]. Due to the biological peculiarities of patients, the actual values differed significantly from the assumptions about liver tissue made in the current model.

To quantify the predictive accuracy of the proposed model, the PINN outputs were compared against the numerical solution generated by the explicit FTCS FDM. The distributions of comparative spatial temperatures and the corresponding absolute error map at $t$ = 600 s are presented in Figure 12.

The global relative $L_2$ error norm for the PINN prediction was evaluated to 1.9419 $\times$ $10^{-2}$ and 4.2079 $\times$ $10^{-2}$ (approximately 1.9% and 4.2% ). The absolute error map indicates that the primary source of this deviation was localized strictly at the geometric center ($x$ = 0, $y$ = 0), corresponding to the focal point of the Gaussian heat source. At this central node, the FDM baseline calculated a maximum temperature of 100.22 °C in both cases, whereas the PINN converged to a peak of 101.51 °C without vascular perfusion and 101.25 °C in the presence of large vessel.

This discrepancy arose from the fundamental differences in how the two methods handled extreme gradients. The explicit FDM advanced through discretized time steps, linearly accumulating thermal energy at the singularity point without the constraints of phase-change mechanics (such as tissue vaporization, which typically caps physical temperatures near 100 °C ). In contrast, PINN acted as a global continuous function approximator. The loss function inherently penalized extreme local gradients to maintain overall continuity and smoothness across the entire 2D + $t$ domain. Consequently, the neural network slightly dampened the non-physical mathematical peak at the singularity. Despite this central deviation, the temperature profiles at the critical 50 °C ablation margin exhibited complete topological agreement, thus confirming that PINN accurately resolved the clinically relevant boundary of the necrosis zone. As seen in Figure 13, the respective maximum and mean absolute errors in PINN predictions are 2.46 °C and 0.76 °C in the absence of heat sink effect against 1.19 °C and 9.32 °C.

One of the main goals of replacing classical numerical solvers with PINN is to reduce inference time, thereby enabling real-time clinical evaluation.

Traditional numerical solvers, such as FEM and FDM, require significant computational time ranging from a few minutes to almost half an hour, depending on the mesh density and dimensionality. Even highly optimized GPU-accelerated methods, such as the Lattice Boltzmann Method, and parallel computing methods, could significantly lead to a longer implementation time to resolve the thermal field [30]. In contrast, the proposed PINN Framework, once trained, reduced the computational load to milliseconds ($<$0.1 seconds), which corresponds to an acceleration factor of approximately $10^{3}$ to $10^{5}$ compared to standard mesh-based solvers.

The total computation time for the FDM solver is directly proportional to the duration of physical time and grid resolution. Conversely, the PINN methodology is based on an offline training phase and an online inference phase. As the trained network effectively stores the parametrized partial differential solution equation, PINN evaluates the spatiotemporal coordinates simultaneously and the inference time to generate a 10 000-point grid for any arbitrary time t is reduced to fractions of a second.

Thanks to the decoupling of physical time progression from computational execution, the PINN framework eliminates sequential time-stepping constraints and provides a rapid prediction mechanism suitable for integration with advanced decision-support systems.

A PINN framework was presented in this study to simulate the spatiotemporal heat transfer dynamics during hepatic RFA. By embedding the two-dimensional Pennes bioheat equation directly into the loss function, the neural network effectively captured the nonlinear interactions between the localized electromagnetic heat source and the convective heat sink effect induced by blood perfusion.

A comparative analysis against a classical explicit FDM yielded a relative $L_2$ error norm of 1.9%. The results demonstrated that while the continuous functional approximation of PINN slightly dampened extreme temperature gradients at the heat source singularity, it accurately resolved the critical 50 °C isothermal boundary required for mapping the coagulation necrosis zone.

The computational evaluation highlighted a structural shift in simulation efficiency. The trained PINN model bypassed the stability constraints that limited traditionally sequential time-stepping algorithms. This decoupling permitted predictions of near-instantaneous temperature field at any arbitrary spatiotemporal coordinate.

Future research should address the physical limitations of the current linear governing equation by incorporating temperature-dependent thermophysical properties and phase-change mechanics to account for tissue vaporization above 100 °C. Furthermore, expanding the computational domain to three dimensions and integrating patient-specific anatomical data will increase the clinical fidelity of the model. Ultimately, the rapid predictive capabilities and physical consistency of PINNs position them as a foundational component for advanced decision-making and simulation technologies in the planning of thermal ablation and broader heat transfer applications.