Boundary-only mechanobiological modeling of thermoelastic tissue deformation during nanoparticle-mediated hyperthermia

Abstract

Localized hyperthermia not only causes thermal injury to the tissues but also causes mechanical deformation of the tissues due to expansion of the tissues resulting from temperature changes, stress creation, and strain concentration. In this paper, an exclusive mechanobiological model for nanoparticle-enhanced hyperthermia is developed by modeling transient bioheat transfer, Arrhenius thermal damage, and thermoelastic responses of the tissues due to temperature change. It combines the transient Pennes bioheat transfer model, spatially heterogeneous nanoparticle-specific absorption rate-based heat sources, Arrhenius thermal damage kinetics, and linear isotropic thermoelasticity based on temperature-based eigenstrain. To circumvent the requirement for domain discretization to incorporate internal heat source, perfusion, and transient heat storage in the bioheat transfer model, a Dual Reciprocity Boundary Element Model is developed. Temperature histories are then used to assess the level of accumulated thermal damage and provide the inputs for thermoelastic boundary integral formulation, where displacements, strains, and stresses will be obtained. Through the correlation of nanoparticle-generated internal heat into mechanical deformations, the proposed methodology establishes the connection between thermal therapy, tissue reaction, and elastography. Numerical validation through manufactured solutions, comparison to the analytical steady state solution, and overall energy balance analysis reveals the first-order convergence and stability of transient simulation along with overall energy conservation up to 2.7% in the transient energy balance test. The simulations show that nanoparticle heating creates detectable thermoelastic displacements, which could be the early signs of treatment.

Graphical Abstract

Introduction

The importance of cancer therapies utilizing thermal ablation and hyperthermia has increased due to their capability to enable energy delivery into malignant tumors without major surgery. This is because radiofrequency ablation, microwave ablation, laser interstitial thermal therapy, focused ultrasound, and nanoparticle-mediated hyperthermia all use the same key principle in their therapeutic approach: an appropriate amount of heat needs to be delivered to the tumor without damaging healthy tissues.^1,2 This requirement makes treatment planning fundamentally a multiphysics problem rather than a purely thermal one. Temperature elevation governs not only heat diffusion, but also thermal dose, protein denaturation, vascular response, tissue stiffness variation, and mechanically detectable deformation. From a thermodynamic and oncological perspective, the ability to connect deposited energy to thermal injury and tissue-scale response is therefore essential for predictive treatment design.^3,4

Classical bioheat transfer models have provided a useful foundation for thermal therapy simulation, but modern biological heat-transfer problems increasingly require formulations that can account for memory effects, non-Fourier conduction, tissue heterogeneity, porous microstructure, and nonlocal interactions. Recent advances in Moore–Gibson–Thompson thermodiffusion, memory-dependent dual-phase-lag bioheat transfer, and coupled magneto-bioheat modeling have shown that generalized heat-transfer theories may be important when classical Fourier diffusion is insufficient to represent thermal relaxation and nonlocal transport in complex biological media.^5–7 These developments are particularly relevant to tumor hyperthermia, where thermal gradients can be steep, material properties can vary spatially, and the response of living tissue may depend on both current and prior thermal exposure.

Boundary-based and generalized computational approaches have also been increasingly used to address bioheat and biothermomechanical problems in complex tissues. Previous boundary element and hybrid numerical models have been developed for anisotropic laser-induced hyperthermia, nonlinear biomechanics of functionally graded biological soft tissues, fractional dual-phase-lag viscoelastic models, anisotropic soft-tissue biomechanics, and time-fractional bioheat transfer in graded tissue media.^8–12 More recent extensions have introduced boundary element frameworks for fractional nanoparticle-assisted photothermal transport, multi-term fractional heat conduction in heterogeneous media, size-dependent micropolar porous thermoelasticity, and generalized thermoelastic sensitivity analysis.^13–16 These studies demonstrate that boundary-integral methods are well suited to problems where high accuracy at boundaries, reduced dimensionality, and treatment of coupled thermal–mechanical fields are important. However, a fully integrated boundary-dominant formulation that links nanoparticle SAR heating, transient perfused bioheat transfer, Arrhenius damage, and thermoelastic deformation remains insufficiently developed.

Magnetic nanoparticles provide a particularly attractive route for localized hyperthermia because they can convert externally applied electromagnetic energy into volumetric heat within nanoparticle-enriched tumor regions. ¹⁷ The classical Rosensweig theory describes the heating of magnetic fluids subjected to alternating magnetic fields and provides a physical basis for representing nanoparticle heating through an effective SAR source term. ¹⁸ In computational oncology, finite element models and imaging-guided simulations have been used to predict magnetic hyperthermia in brain tumors, liver tumors, and magnetic-particle-imaging-guided treatment scenarios.^19–21 These studies show the clinical relevance of nanoparticle heating, but they are generally domain-discretization based and often focus primarily on temperature prediction and thermal dose. A model that also predicts deformation, strain, and stress from the same thermal field can provide a more complete description of tissue response during treatment.

The Pennes bioheat equation remains the most widely used starting point for perfused tissue heat-transfer modeling because it provides a tractable continuum representation of conduction, perfusion-mediated heat exchange, and volumetric heat generation. ²² Extensions of tissue heat-transfer theory have also emphasized the role of vascular microstructure, thermal relaxation, and coupled biothermomechanics in determining local temperature fields and tissue response.^23–25 Once the temperature field is known, tissue injury is often evaluated using thermal damage models based on temperature history. Theoretical and experimental studies of thermal injury have established that irreversible damage depends not only on peak temperature, but also on exposure duration.^26–28 Related imaging and biochemical studies have further shown that heat-shock and denaturation responses can be connected to measurable thermal histories. ²⁹ Therefore, any model intended for treatment planning should compute damage as a time-integrated local field rather than as a simple threshold based only on instantaneous temperature.

In parallel, elastography has become an important technique for assessing tissue mechanical properties and deformation. Ultrasound and related elastographic methods provide noninvasive access to displacement and strain information, making them potentially valuable for monitoring thermal therapy.^30–32 First, in nanoparticle-assisted hyperthermia treatment, thermal expansion arises owing to the heating effect. When this expansion is non-uniform or mechanically restrained due to interaction with neighboring tissue, displacements, strains, and stresses occur. These continuum mechanical fields do not by themselves constitute a cell-scale mechanobiological model, but they provide mechanical indicators of the evolving thermal state and can help interpret externally measurable deformation. Thus, a coupled thermal–damage–thermoelastic model can provide a bridge between internal nanoparticle heating and potential elastographic monitoring.

Boundary Element Methods offer a powerful alternative to domain-based discretization for such problems because they reduce the dimensionality of the computational domain and provide accurate boundary and interface representations. Classical boundary integral equation methods have been widely developed for solid and fluid mechanics, fracture mechanics, and thermo-fluid problems.^33–37 The Dual Reciprocity Boundary Element Method extends the applicability of BEM to problems containing domain terms by approximating nonhomogeneous sources through basis functions and corresponding particular solutions. ³⁸ This is crucial for transient bioheat transfer because the governing equation contains volumetric terms associated with transient storage, perfusion, SAR heating, and thermal eigenstrain. Energy-based verification concepts and transient heat-conduction DR-BEM formulations further support the suitability of boundary-integral methods for time-dependent multiphysics transport problems.^39,40 Recent fractional and memory-driven boundary element formulations also indicate that BEM-based approaches can be extended to increasingly general bioheat models. ⁴¹ Earlier applications of BEM to the bioheat equation and transient diffusion–convection–reaction problems provide additional precedent for using boundary-integral formulations in biological heat-transfer analysis.^42,43

Despite these advances, several gaps remain. First, many nanoparticle hyperthermia models focus on temperature alone and do not explicitly propagate the thermal solution into damage and mechanical deformation. Second, most thermoelastic models of hyperthermia use analytical solutions without consistent coupling between the temperature and mechanical field. Third, the predominant computational tools are still the finite element and finite difference method, which both need a volumetric mesh and thus could be costly for parameter studies with localized sources. Fourth, existing boundary-based bioheat formulations often require clearer treatment of internal source terms, interface continuity, thermal eigenstrain, and validation against analytical or finite element benchmarks.

In this paper, a coupled boundary-dominant computational framework is developed for nanoparticle-mediated hyperthermia in perfused biological tissue. The proposed model combines a Pennes-type transient bioheat equation with a spatially distributed nanoparticle SAR source, Arrhenius thermal damage kinetics, and linear isotropic thermoelasticity driven by temperature-dependent eigenstrain. The nanoparticle heat source is represented as a smooth radially decaying SAR distribution, while the accumulated thermal damage is evaluated pointwise from the local temperature history. The thermoelastic response is then computed from the temperature-induced eigenstrain field, allowing prediction of displacement, strain, and stress as continuum-scale mechanical indicators of heating.

The numerical formulation is based on the Dual Reciprocity Boundary Element Method. At each time step, the transient bioheat equation is reformulated as a Poisson-type boundary integral problem. Domain terms arising from transient storage, perfusion, and SAR heating are approximated using Dual Reciprocity basis functions and converted into equivalent boundary integrals through analytical particular solutions. The resulting temperature field is then used to update the Arrhenius damage integral and to drive the thermoelastic boundary element formulation. Although the method avoids volumetric meshing, it is described here as boundary-dominant rather than strictly boundary-exclusive, because a limited set of interior collocation points is used for Dual Reciprocity approximation, damage evaluation, and field reconstruction.

The main contribution of this work is the establishment of a consistent thermal-to-mechanical computational pathway for nanoparticle-mediated hyperthermia. The framework links localized SAR heating to transient temperature evolution, temperature history to cumulative thermal damage, and temperature gradients to thermoelastic displacement, strain, and stress. This provides a physically interpretable continuum model for studying how internal heating may produce mechanically observable signals at tissue boundaries. The formulation is verified using manufactured solutions, analytical heat-conduction comparisons, time-step and mesh refinement, global energy-balance assessment, and an independent finite element benchmark for selected thermoelastic quantities. The resulting model provides a foundation for future extensions involving fractional or non-Fourier heat transfer, temperature-dependent tissue properties, viscoelasticity, stress-dependent perfusion, experimentally measured nanoparticle distributions, and elastography-based validation.

Governing physical framework

The goal of this section is to provide a consistent physical framework that connects three phenomena that are all too often considered separately: nanoparticle-induced energy deposition, non-uniform heat transport in a living medium, and deformations caused by thermal expansion. In reality, these phenomena cannot be separated. A localized heating event necessarily affects temperature, causes biological damage, and produces stress and motion in the surrounding medium at the same time. A physical model that does not take any of these phenomena into account is physically and biologically incomplete.

Rather than combining disparate equations, the current framework is based on a chain of cause-and-effect relationships. The energy imparted by magnetic/optical excitation of nanoparticles is represented as a volumetric heat source. This heat source is used to describe the transient evolution of temperature in perfused tissue, which is used to determine irreversible tissue damage using kinetic damage laws, and spatial temperature gradients are used to determine eigenstrains that cause propagating tissue deformation, represented on the tissue surface by elastographic imaging modalities.

Two design principles are used in the formulation. First, each equation must relate to a quantity that is measurable or inferable: temperature, heat flux, damage, stress, or displacement. Second, each coupling must relate to a physical process, not a numerical convenience. Heating couples with mechanics through thermal expansion; mechanics does not feedback to affect heat transport in the current formulation because of the time scale separation observed in soft tissue.

The domain Ω describes a tissue region with a nanoparticle-enriched tumor subdomain embedded in healthy host tissue. The boundary Γ can be subject to temperature, convective, and adiabatic boundary conditions. The internal heat source is provided solely by the nanoparticles.

Within this setting, four coupled components are introduced in the following subsections: (i) a nanoparticle-based volumetric heat source, (ii) a transient bioheat equation with perfusion, (iii) a kinetic thermal damage model, and (iv) a thermoelastic formulation linking temperature to deformation.

Together, these components define a unified magneto–thermo–mechanical description of localized thermal therapy, designed not only to predict temperature, but to translate heating into biologically and mechanically meaningful outcomes.

Geometrical configuration and tissue domains

The physical setting considered in this work is a tumor embedded within a surrounding region of healthy biological tissue, representative of a localized malignancy treated by nanoparticle-mediated magnetothermal therapy. The computational domain is idealized as a circular (in two dimensions) or spherical (in three dimensions) tissue region of outer radius $R_{\text{ext}}$ , containing a centrally located tumor of radius $R_t$ . This geometry reflects key characteristics of clinically relevant configurations while allowing a clear and rigorous investigation of coupled thermal and mechanical effects.

The tissue body is decomposed into two subdomains: the tumor region ${\mathrm{\Omega }}_t$ and the surrounding healthy tissue ${\mathrm{\Omega }}_h$ , such that

$\mathrm{\Omega }={\mathrm{\Omega }}_t\cup {\mathrm{\Omega }}_h,{\mathrm{\Omega }}_t\cap {\mathrm{\Omega }}_h=\varnothing ,$

(1)

with an internal interface ${\mathrm{\Gamma }}_t=\partial {\mathrm{\Omega }}_t$ and an external boundary $\mathrm{\Gamma }=\partial \mathrm{\Omega }$ . The tumor tissue interface, denoted by the set ${\mathrm{\Gamma }}_t$ is taken to be perfectly bonded, which means the temperature, heat flux, displacement, and traction are all continuous across the interface.

More explicitly, at the tumor–host interface Γ_th, the thermal and mechanical continuity conditions are imposed as

At the tumor–host interface Γ_th​:

$T_t=T_h,k_t\frac{\partial T_t}{\partial n}=k_h\frac{\partial T_h}{\partial n},$

(2)

${\mathbf{u}}_t={\mathbf{u}}_h,{\mathbf{t}}_t+{\mathbf{t}}_h=\mathbf{0},$

(3)

For thermal boundary:

$-k\frac{\partial T}{\partial n}=h\left(T-T_{\infty }\right).$

(4)

At the external boundary:

$\mathbf{t}=\mathit{\sigma }\mathbf{n}=\mathbf{0}.$

(5)

where the subscripts $t$ and $h$ denote tumor and host tissue, respectively. These conditions enforce continuity of temperature, continuity of normal heat flux, continuity of displacement, and balance of traction across the perfectly bonded interface.

This is because there is no mechanical separation between the tumor and the host tissue, and the use of any kind of artificial thermal or mechanical barriers is avoided.

Due to the presence of a centrally located heat source and a radially symmetric concentration distribution of nanoparticles in hyperthermia treatment protocols, axisymmetric conditions are employed. In addition, the spatial coordinates are represented in polar form $\left(r,\theta \right)$ for a 2D problem or spherical coordinates $\left(r,\theta ,\varphi \right)$ for a 3D problem, with all field variables being functions only of the radial coordinate r and time t.

The boundary $\mathrm{\Gamma }$ describes the boundary between the modeled region of tissue and the surrounding physiological environment. In terms of heat transfer, this boundary can have prescribed temperature (Dirichlet), convective (Robin), and/or mixed conditions, which describe the cooling effects of the vascular network and the interaction with the surrounding tissue layers. In terms of the mechanical boundary conditions, the boundary $\mathrm{\Gamma }$ is assumed to have traction-free conditions, unless otherwise specified, to describe a region of soft tissue that is not subject to any rigid anatomical structure. This boundary condition allows for thermally induced expansion and displacement of the tissue, which is important for elastographic interpretation.

This geometric setup provides a clean and physically consistent setting in which volumetric nanoparticle heating is localized in ${\mathrm{\Omega }}_t$ , thermal energy transport occurs across $\mathrm{\Omega }$ , and thermoelastic deformation is imparted to the external boundary $\mathrm{\Gamma }$ . It is a minimal but realistic setting for exploring the impact of localized magnetothermal energy delivery on temperature increase, tissue injury, and mechanically detectable motion on the boundary.

Bioheat transfer with nanoparticle heating

The heat transfer through the tissue region is modeled using an extended version of the Pennes bioheat equation, where the volume heating effect of magnetic nanoparticles is included explicitly. The notation used is $T\left(\mathbf{x},t\right)$ as the temperature at location $\mathbf{x}$ inside Ω and at time $t$ , $\rho$ being the density, $c$ the specific heat, and k the thermal conductivity of the tissue.

The transient thermal balance in the tissue is written as

$\rho c\frac{\partial T}{\partial t}=\nabla ·\left(k\nabla T\right)-{\rho }_b{c}_b{\omega }_b\left(T-T_b\right)+Q_{\text{SAR}}\left(\mathbf{x},t\right),\mathbf{x}\in \mathrm{\Omega },$

(6)

where the notation ${\rho }_b$ , $c_b$ , and ${\omega }_b$ are used to represent the density of the blood, the specific heat of the blood, and the volumetric perfusion rate, respectively, while $T_b$ is the arterial blood temperature. The perfusion term describes the convective heat transfer process between the tissue and the flowing blood. This term acts as a distributed sink and can have different values for tumor and healthy tissues. Perfusion in tumor tissue is usually lower than in healthy tissue.

The source term $Q_{\text{SAR}}\left(\mathbf{x},t\right)$ models energy deposition by magnetic nanoparticles subjected to an external alternating magnetic field. It represents the local specific absorption rate (SAR) and is expressed as

$Q_{\text{SAR}}\left(\mathbf{x},t\right)=A\hspace{0.17em}S\left(t\right)\hspace{0.17em}\chi \left(\mathbf{x}\right),$

(7)

where $A$ is the peak volumetric SAR amplitude, $S\left(t\right)$ describes the temporal heating protocol, and $\chi \left(\mathbf{x}\right)$ is the normalized spatial distribution of nanoparticle absorption. Since the present configuration is radially symmetric, the spatial function $\chi \left(\mathbf{x}\right)$ reduces to $\chi \left(r\right)$ , where $r=\Vert \mathbf{x}-{\mathbf{x}}_c\Vert$ , is the distance from the tumor center, and ${\mathbf{x}}_c$ denotes the tumor center. In the baseline simulations, a smooth Gaussian distribution is adopted:

$\chi \left(r\right)=\exp \left[-{\left(\frac{r}{r_0}\right)}^2\right]$

(8)

where $r_0$ is the SAR decay radius. Therefore, the explicit nanoparticle heating source becomes

$Q_{\text{SAR}}\left(r,t\right)=A\hspace{0.17em}S\left(t\right)\exp \left[-{\left(\frac{r}{r_0}\right)}^2\right]$

(9)

This form provides a continuous localized heating profile across the tumor–host interface and avoids an artificial jump in the volumetric heat source. During the heating interval, the temporal factor is taken as $S\left(t\right)=1$ which corresponds to a cycle-averaged representation of the alternating magnetic field. Thus, the rapid magnetic-field oscillations are not resolved explicitly; instead, their thermal effect is incorporated into the effective SAR amplitude A.

The initial condition is prescribed as a uniform physiological temperature,

$T\left(\mathbf{x},0\right)=T_0,\mathbf{x}\in \mathrm{\Omega },$

(10)

where $T_0$ ​is typically taken as the normal core body temperature. At the outer boundary $\mathrm{\Gamma }$ , heat exchange with surrounding tissue and blood flow is represented using prescribed temperature conditions.

$T=T\prime \hspace{1em}\text{on}\hspace{0.17em}{\mathrm{\Gamma }}_T,$

(11)

or convective exchange,

$-k\nabla T·\mathbf{n}=h\left(T-T_{\infty }\right)\hspace{1em}\text{on}\hspace{0.17em}{\mathrm{\Gamma }}_q,$

(12)

where $h$ is the heat transfer coefficient, $T_{\infty }$ ​ is the reference temperature of the surrounding and nis the outward direction normal to the boundary. All these terms combined give the physiological cooling effect through which heat can exit the system.

This equation summarizes the physics behind hyperthermia using nanoparticles; this includes the volumetric heating effect in the tumor, diffusion of heat to the surrounding tissues, and continuous extraction of energy through perfusion. It is the thermal foundation of the coupled model, from which both the biological effects and the mechanical response are subsequently derived.

Arrhenius thermal damage model

Elevation of temperature through nanoparticle-assisted hyperthermia might lead to temporary thermal stress at the beginning stage in biological tissue. In cases where the temperature is sustained high enough for a long period of time, irreversible biochemical and structural changes can take place, such as protein denaturation and inactivation of enzymes. Therefore, thermal injury cannot be determined from the instantaneous temperature alone. It must be evaluated from the complete local temperature history.

In the present formulation, cumulative thermal injury is quantified using the Arrhenius damage model. At any material point $\mathbf{x}$ and time $t$ , the dimensionless thermal damage integral is defined as

$\mathrm{\Omega }\left(\mathbf{x},t\right)={\int }_0^{\mathrm{t}} A_0\hspace{0.17em}\exp \left[-\frac{E_{\mathrm{a}}}{R\hspace{0.17em}T\left(\mathbf{x},\tau \right)}\right]d\tau .$

(13)

where $\mathrm{\Omega }\left(\mathbf{x},t\right)$ is the accumulated thermal damage, $A_0$ is the frequency factor, $E_{\mathrm{a}}$ is the activation energy for thermal denaturation, R is the universal gas constant, and $T\left(\mathbf{x},\tau \right)$ is the absolute temperature at position $\mathbf{x}$ and previous time $\tau$ . This expression reflects the exponential sensitivity of biological damage to temperature. Small temperature increases may have limited effect, whereas modest increases in the hyperthermic or ablative range can substantially accelerate tissue damage accumulation.

The local survival fraction is related to the Arrhenius damage integral by

$S_{\mathrm{v}}\left(\mathbf{x},t\right)=\exp \left[-\mathrm{\Omega }\left(\mathbf{x},t\right)\right].$

(14)

where $S_{\mathrm{v}}\left(\mathbf{x},t\right)$ denotes the fraction of surviving tissue or cells at point x and time t. Accordingly, the local fraction of damaged tissue is

$D\left(\mathbf{x},t\right)=1-\exp \left[-\mathrm{\Omega }\left(\mathbf{x},t\right)\right].$

(15)

In this study, $\mathrm{\Omega }=1$ is adopted as the threshold for significant irreversible thermal injury. At this threshold,

$D=1-\exp \left(-1\right)\approx 0.632$

(16)

which corresponds to approximately 63% thermal damage or cell death. Thus, regions satisfying $\mathrm{\Omega }\left(\mathbf{x},t\right)\ge 1$ are interpreted as significantly damaged, whereas regions with $\mathrm{\Omega }\left(\mathbf{x},t\right)<1$ are interpreted as sub-lethally heated or reversibly affected, depending on the local exposure history.

Numerically, the Arrhenius damage integral is evaluated pointwise at each interior collocation point used for field reconstruction. It is not computed from a tumor-volume averaged temperature. At a collocation point xᵢ, the discrete damage update is written as

${\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{n}}={\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{n}-1}+\mathrm{\Delta }t\hspace{0.17em}{A}_0\hspace{0.17em}\exp \left[-\frac{E_{\mathrm{a}}}{\left(R{T}_{\mathrm{i}}^{\mathrm{n}}\right)}\right].$

(17)

where ${\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{n}}$ is the damage at collocation point xᵢ and time level $t^{\mathrm{n}},T_{\mathrm{i}}^{\mathrm{n}}$ is the local absolute temperature obtained from the DR-BEM thermal solution, and $\mathrm{\Delta }t$ is the time step. A trapezoidal update may also be used for improved temporal integration accuracy:

${\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{n}}={\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{n}-1}+\left(\frac{\mathrm{\Delta }t}{2}\right)\left\{A_0\hspace{0.17em}\exp \left[-\frac{E_{\mathrm{a}}}{\left(R{T}_{\mathrm{i}}^{\mathrm{n}-1}\right)}\right]+A_0\hspace{0.17em}\exp \left[-\frac{E_{\mathrm{a}}}{\left(R{T}_{\mathrm{i}}^{\mathrm{n}}\right)}\right]\right\}.$

(18)

This pointwise implementation is essential because nanoparticle-mediated heating is spatially non-uniform. Two points within the same tumor may experience different thermal histories and therefore accumulate different damage levels even when they belong to the same anatomical region. The resulting spatial field $\mathrm{\Omega }\left(\mathbf{x},t\right)$ provides a direct map of the predicted thermal injury zone.

The damage field reported in the results is evaluated from the local temperature history at each collocation point. It is not computed from a prescribed radial damage profile or from a tumor-averaged temperature. For each spatial point, the full transient temperature sequence $T\left(\mathbf{x},t\right)$ is inserted into the Arrhenius integral using the parameters $A_0\hspace{0.17em}\text{and}\hspace{0.17em}{E}_{\mathrm{a}}$ listed in Table 1. This ensures that spatial variations in $\mathrm{\Omega }$ arise directly from the computed thermal field and the exposure duration.

Table 1.

Arrhenius thermal injury parameters and damage criteria adopted in the simulations.

Parameter	Symbol	Value used	Unit	Role in the model
Frequency factor	$A_0$	$3.1\times {10}^{98}$	s−1	Controls the pre-exponential rate of thermal injury accumulation.
Activation energy	$E_{\mathrm{a}}$	$6.28\times {10}^5$	J mol−1	Defines the temperature sensitivity of tissue denaturation.
Universal gas constant	$R$	8.314	J mol-1 K-1	Physical constant used in the Arrhenius exponential term.
Local absolute temperature	$T\left(\mathbf{x},\tau \right)$	Computed from DR-BEM	K	Local transient temperature history used to evaluate thermal damage.
Arrhenius damage integral	$\mathrm{\Omega }\left(\mathbf{x},t\right)$	Computed pointwise	—	Cumulative measure of thermal injury at each interior collocation point.
Survival fraction	$S_{\mathrm{v}}\left(\mathbf{x},t\right)$	$\exp \left[-\mathrm{\Omega }\left(\mathbf{x},t\right)\right]$	—	Fraction of tissue or cells remaining viable after thermal exposure.
Damage fraction	$D\left(\mathbf{x},t\right)$	$1-\exp \left[-\mathrm{\Omega }\left(\mathbf{x},t\right)\right]$	—	Fraction of thermally damaged tissue.
Irreversible injury threshold	$\mathrm{\Omega }\text{crit}$	1	—	Criterion for significant irreversible thermal injury.
Damage fraction at threshold	$D\text{crit}$	$1-\exp \left(-1\right)\approx 0.632$	—	Corresponds to approximately 63% thermal damage or cell death.
Time step	$\mathrm{\Delta }t$	$0.5$	s	Time increment used to update the Arrhenius damage integral.

The Arrhenius parameters used in the simulations are summarized in Table 1. These parameters are tissue-dependent and should be selected from experimentally reported thermal damage data for the target tissue type. In the present proof-of-concept simulations, representative soft-tissue values are adopted to demonstrate the coupled thermal–damage–mechanical framework. The sensitivity of Ω to A₀ and E_a should be recognized, since the exponential form of the Arrhenius law can amplify parameter uncertainty at elevated temperatures.

Although the present work does not perform direct experimental validation of thermal damage against measured histological injury, the use of the Arrhenius model provides a standard and physically interpretable link between transient temperature exposure and irreversible thermal injury. The damage field is therefore used as a mechanobiological indicator of treatment progression and as a biological context for interpreting the computed thermoelastic displacement, strain, and stress fields.

Thermoelastic tissue response

Thermal loading of biological tissue is intrinsically mechanical. Upon heating, tissues expand volumetrically. When this expansion is heterogeneous, as with nanoparticle-based heating, stress and macroscopic deformation result. These thermoelastic phenomena are not ancillary; they are the cause of the surface displacements measured by elastography during thermal therapy. To model this phenomenon, the biological tissues are assumed to behave like a linear, isotropic thermoelastic material undergoing small deformations.

The infinitesimal strain tensor can be written as

$\mathit{\epsilon }=\frac{1}{2}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\mathrm{T}}\right).$

(19)

Thermal expansion enters the kinematics through an eigenstrain

${\mathit{\epsilon }}^{\theta }=\alpha \left(T-T_0\right)\mathbf{I},$

(20)

where $\alpha$ is the coefficient of linear thermal expansion, $T_0$ a reference temperature, and $\mathbf{I}$ the identity tensor. The elastic strain driving stress is therefore

${\mathit{\epsilon }}^e=\mathit{\epsilon }-{\mathit{\epsilon }}^{\theta }.$

(21)

For an isotropic medium characterized by Young’s modulus $E$ and Poisson’s ratio $\nu$ , the Cauchy stress is given by

$\mathit{\sigma }=\mathbf{C}:{\mathit{\epsilon }}^e=\mathbf{C}:\left(\mathit{\epsilon }-\alpha \left(T-T_0\right)\mathbf{I}\right),$

(22)

where $\mathbf{C}$ is the fourth-order elasticity tensor, and equilibrium is assumed in the quasi-static regime.

$\nabla ·\mathit{\sigma }=\mathbf{0},\mathbf{x}\in \mathrm{\Omega }.$

(23)

The outer boundary $\mathrm{\Gamma }$ is taken to be traction-free unless specified otherwise.

$\mathit{\sigma }\mathbf{n}=\mathbf{0}\hspace{0.17em}\text{on}\hspace{0.17em}\mathrm{\Gamma },$

(24)

This allows the tissue to expand and deform due to heating. The continuity of displacement and force is imposed on the interface between the tumor and tissue, denoted by ${\mathrm{\Gamma }}_t$ , This considers the mechanical binding.

In this manner, it is possible to relate the thermal and mechanical fields by the eigenstrain term. The temperature distribution, solved using the bioheat equation, is used as a body force in the mechanics problem. Regions of high rates of heating will lead to differential expansion, which is constrained by the surrounding tissues. This will lead to stress and consequent deformation away from the heated region. Notably, this deformation is not restricted to the heated region but will propagate towards the boundary, where it can be measured using elastography.

From a mechanobiological perspective, the resulting thermoelastic strain should not be regarded merely as a secondary mechanical byproduct of heating. Instead, it represents a mechanically encoded signature of the evolving thermal and biological state of the tissue. Localized temperature gradients generate spatially heterogeneous expansion, and the surrounding tissue converts this expansion into measurable strain and stress fields. These mechanical fields may carry information about internal energy deposition, tissue constraint, perfusion-mediated cooling, and the onset of thermal damage. Therefore, thermoelastic deformation provides a direct mechanobiological link between nanoparticle-induced heating and externally observable tissue-scale mechanical behavior.

Within the present framework, the contribution of thermoelasticity is twofold. On the one hand, it models the stress and strain values that might potentially influence perfusion, microstructure, and the post-processing tissue response. On the other hand, it provides a deterministic mapping between the internal heating and the resulting surface motion. By including this thermoelastic model in the bioheat formulation, the deformation process will no longer be a qualitative outcome but a predictive field variable.

Mechanobiological interpretation of thermoelastic fields

The thermoelastic fields computed in this work have a direct mechanobiological interpretation. Temperature-induced eigenstrain represents the local tendency of tissue to expand in response to heating, whereas the resulting stress field reflects the mechanical constraint imposed by neighboring tissue regions. In nanoparticle-mediated hyperthermia, this distinction is essential. Localized heating can occur due to high absorption of energy within the tumor region; however, expansion of the heated area is opposed by the normal tissues in the vicinity. This interaction generates spatially distributed strain and stress fields that encode the mechanical consequences of the thermal treatment.

From a mechanobiological viewpoint, these fields are not merely numerical outputs. Displacement, strain, and stress are mechanical cues that may influence the tumor microenvironment. These local compressive and tensile forces can affect vascular patency, modify interstitial pressure, modify the efficiency of perfusion, and modify heat removal due to blood flow. For the mechanically altered region of thermal damage, the mechanical behavior may change depending on changes in tissue stiffness and fluid permeability due to tissue damage. The thermoelastic response, therefore, acts as a mechanical link to thermal effects on biological systems.

The temperature-induced strain field is especially important because it appears immediately as temperature gradients develop. Unlike cumulative thermal damage, which depends on the full exposure history, thermoelastic strain is an instantaneous mechanical signature of the evolving temperature field. This makes it potentially useful as an early indicator of energy localization and treatment progression. If measured by ultrasound elastography, optical coherence elastography, magnetic resonance elastography, or related techniques, the displacement field may provide indirect information about the internal heating pattern before irreversible tissue injury is fully established.

Tumorous tissue has an abnormal vascular structure as well as varied stiffness, which makes it particularly sensitive to mechanical loading due to these structural and material characteristics. Thermal stresses can be among the factors involved in causing vascular compression, perfusion modification, or interstitial transport. Although the present formulation treats perfusion as a prescribed thermal parameter, the computed stress field provides a natural basis for future coupling between mechanics and blood flow. Such coupling would allow the model to represent stress-dependent perfusion, vascular collapse at high temperature, and damage-induced modification of tissue transport properties.

The present formulation does not explicitly solve cellular mechanotransduction pathways. Although the mechanism does not include the role of mechanosensitive ion channels, cytoskeletal reorganization, extracellular matrix restructuring, or biochemical signal transmission on the cellular level, it offers the continuum-level mechanics which form the input for these mechanisms. The displacement, strain, and stress fields predicted by this approach constitute the mechanical environment during hyperthermia. These quantities can therefore serve as a bridge between macroscopic treatment planning and future models of cellular or tissue-level mechanobiological response.

In this sense, the proposed model connects bioheat transfer and mechanobiology. Temperature is not treated as the final endpoint of the analysis, but as the driver of mechanical cues that can be observed, quantified, and potentially used for treatment monitoring. The framework translates nanoparticle-mediated heating into tissue deformation, strain localization, and stress generation. This thermal-to-mechanical pathway provides a mechanobiological basis for elastographic feedback and supports the broader interpretation of hyperthermia as both a thermal and mechano-therapeutic intervention.

Boundary Element–Dual Reciprocity formulation

In many bioheat models, the main limitation comes from the numerical treatment rather than the underlying physics. In many bioheat models, the main limitation comes from the numerical treatment rather than the underlying physics.

In many bioheat models, the main limitation comes from the numerical treatment rather than the underlying physics.

In many bioheat models, the main limitation comes from the numerical treatment rather than the underlying physics. The actual use of the volume method will become vulnerable when the heat source is local, the volume is large, and biological interfaces are present. After years of wrestling with mesh distortion near tumors and re-meshing whenever a geometry changes, it becomes clear that the problem itself is boundary-dominated: what clinicians care about is what happens at interfaces and at the surface.

The Boundary Element Method (BEM) is therefore a natural vehicle for this class of problems. Its classical limitation—the presence of domain terms in transient and perfused bioheat equations—is overcome here through the Dual Reciprocity (DR) technique, which converts volumetric effects into equivalent boundary contributions. The result is a genuinely boundary-only formulation for a fully coupled thermal–mechanical system.

Boundary integral representation of the bioheat equation

The transient temperature field is governed by the Pennes bioheat equation with a spatially distributed nanoparticle heat source. In the present numerical implementation, the computational domain is treated as a two-dimensional circular tissue section under radial symmetry. Therefore, the boundary integral formulation is constructed using the two-dimensional fundamental solution of the Laplace operator. The corresponding three-dimensional kernel is also stated for completeness and to clarify the natural extension to spherical geometries; however, all numerical results reported in this work are obtained using the two-dimensional logarithmic kernel.

At the discrete time level $t^{\mathrm{n}},$ the transient bioheat equation is advanced using the backward Euler scheme. The governing equation is written in the following Poisson-type form:

${\nabla }^2{T}^{\mathrm{n}}\left(\mathbf{x}\right)=b^{\mathrm{n}}\left(\mathbf{x}\right).$

(25)

where Tⁿ(x) is the temperature at the current time step and bⁿ(x) contains all nonhomogeneous domain contributions, namely transient storage, blood perfusion, and nanoparticle-induced SAR heating. The source function is defined as

$b^{\mathrm{n}}\left(\mathbf{x}\right)=\frac{\left(\frac{\rho c}{k}\right)\left[T^{\mathrm{n}}\left(\mathbf{x}\right)-T^{\mathrm{n}-1}\left(\mathbf{x}\right)\right]}{\mathrm{\Delta }t}+\left(\frac{{\rho }_{\beta }{c}_{\beta }{\omega }_{\beta }}{k}\right)\left[T^{\mathrm{n}}\left(\mathbf{x}\right)-T_{\beta }\right]-\left(\frac{1}{k}\right)Q_{\text{SAR}}\left(\mathbf{x},t^{\mathrm{n}}\right)$

(26)

where $\rho$ , $c$ , and $k$ stand for the density of the tissue, the specific heat, and the thermal conductivity of the tissue, respectively. ${\rho }_{\beta }$ , $c_{\beta }$ , and ${\omega }_{\beta }$ stand for the density of the blood, the specific heat of the blood, and the perfusion rate, respectively. In addition, $T_{\beta }$ is the arterial blood temperature, $\mathrm{\Delta }t$ is the time step, and $Q_{\text{SAR}}$ is the nanoparticle-mediated volumetric heating source.

For the two-dimensional circular domain used in the present simulations, the Green’s function of the Laplace operator is given by

$G\left(\mathbf{x},\mathit{\xi }\right)=-\left(\frac{1}{2\pi }\right)\ln \hspace{0.17em}r,r=\Vert \mathbf{x}-\xi \Vert .$

(27)

where $\mathit{\xi }$ is the source point, $\mathbf{x}$ is the field point, and $r$ is the Euclidean distance between them. The outward normal derivative of the Green’s function is

$\frac{\partial \mathrm{G}\left(\mathbf{x},\mathrm{\xi }\right)}{\partial \mathrm{n}}=-\frac{\left(\frac{1}{2\mathrm{\pi }}\right)\left[\left(\mathbf{x}-\mathrm{\xi }\right)·\mathbf{n}\right]}{{\mathrm{r}}^2}.$

(28)

where $\mathbf{n}$ is the outward unit normal vector at the boundary point.

Applying Green’s second identity to the Poisson form of the bioheat equation gives the boundary integral representation

$c\left(\xi \right)T^{\mathrm{n}}\left(\xi \right)+{\int }_{\mathrm{\Gamma }}{T}^{\mathrm{n}}\left(\mathbf{x}\right)\left[\frac{\partial G\left(\mathbf{x},\xi \right)}{\partial n}\right]d\mathrm{\Gamma }\left(\mathbf{x}\right)$

$={\int }_{\mathrm{\Gamma }}G\left(\mathbf{x},\xi \right)q^{\mathrm{n}}\left(\mathbf{x}\right)d\mathrm{\Gamma }\left(\mathbf{x}\right)+{\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right)b^{\mathrm{n}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right).$

(29)

where $\mathrm{\Gamma }$ is the boundary of the tissue domain, $\mathrm{\Omega }$ is the tissue domain, and $c\left(\xi \right)$ is the free-term coefficient associated with the position of the source point. The quantity $q^{\mathrm{n}}\left(\mathbf{x}\right)$ denotes the outward normal temperature gradient and is defined as

$q^{\mathrm{n}}\left(\mathbf{x}\right)=\frac{\partial T^{\mathrm{n}}\left(\mathbf{x}\right)}{\partial n}.$

(30)

If the physical heat flux is required, it is obtained from Fourier’s law as

$q^n\left(\mathbf{x}\right)=-k\frac{\partial T^n\left(\mathbf{x}\right)}{\partial n}.$

(31)

The final domain integral in the boundary integral representation arises from transient storage, perfusion, and internal SAR heating. This integral prevents the direct use of a purely classical boundary-only formulation. In the present work, this term is transformed into equivalent boundary contributions using the Dual Reciprocity Method. Thus, the present subsection establishes the Poisson-type boundary integral form, while the following subsection explains the conversion of the remaining domain term into boundary-only contributions.

For completeness, the corresponding three-dimensional Green’s function for a spherical tissue domain is

$G_3\mathrm{D}\left(\mathbf{x}.\xi \right)=\frac{1}{4\pi r},r=\Vert \mathbf{x}-\mathit{\xi }\Vert .$

(32)

and its outward normal derivative is

$\frac{\partial G_3\mathrm{D}\left(\mathbf{x},\xi \right)}{\partial n}=-\frac{\left[\left(\mathbf{x}-\xi \right)·\mathbf{n}\right]}{\left(4\pi r^3\right)}.$

(33)

This three-dimensional kernel is not used in the two-dimensional simulations reported in this work. It is included only to show that the same DR-BEM framework can be extended naturally to spherical or fully three-dimensional patient-specific geometries. Therefore, the present formulation should be interpreted as a two-dimensional boundary-only implementation, with the three-dimensional form retained as a direct mathematical extension.

Dual Reciprocity approximation of domain terms

The boundary integral equation derived in Section 3.1 contains a domain integral associated with the nonhomogeneous term $b^{\mathrm{n}}\left(\mathbf{x}\right).$ This term contains the effects of transient heat storage, blood perfusion, and nanoparticle-induced SAR heating. In a classical Boundary Element Method, such a domain integral would require internal domain discretization, thereby removing the boundary-only advantage of the method. To preserve the boundary-only structure, the present formulation uses the Dual Reciprocity Method to approximate the domain contribution and transform it into equivalent boundary integrals.

After backward Euler time discretization, the right-hand-side function of the Poisson-type bioheat equation is written as

$b^{\mathrm{n}}\left(\mathbf{x}\right)=\frac{\left(\frac{\rho c}{k}\right)\left[T^{\mathrm{n}}\left(\mathbf{x}\right)-T^{\mathrm{n}-1}\left(\mathbf{x}\right)\right]}{\mathrm{\Delta }t}+\left(\frac{{\rho }_{\beta }{c}_{\beta }{\omega }_{\beta }}{k}\right)\left[T^{\mathrm{n}}\left(\mathbf{x}\right)-T_{\beta }\right]-\left(\frac{1}{k}\right)Q_{\text{SAR}}\left(\mathbf{x},t^{\mathrm{n}}\right).$

(34)

where the first term represents transient thermal storage, the second term represents perfusion-mediated heat exchange with blood, and the third term represents the localized heat input produced by magnetic nanoparticles. The objective of the Dual Reciprocity approximation is to represent this nonhomogeneous function by a finite expansion of known basis functions. Thus,

$b^n\left(\mathbf{x}\right)\approx {\displaystyle \sum }_{j=1}^N {\alpha }_j^n{\varphi }_j\left(\mathbf{x}\right).$

(35)

where

$N=N_{\beta }+N_{\mathrm{i}}.$

(36)

Here, $N_{\beta }$ is the number of boundary nodes, $N_{\mathrm{i}}$ is the number of interior collocation points, ${\varphi }_{\mathrm{j}}\left(\mathbf{x}\right)$ are prescribed radial basis or polynomial basis functions, and ${\alpha }_j^n$ are time-dependent Dual Reciprocity coefficients at the current time step.

For each basis function ${\varphi }_j\left(\mathbf{x}\right)$ , an analytical particular solution ${\psi }_{\mathrm{j}}\left(\mathbf{x}\right)$ is introduced such that

${\nabla }^2{\psi }_{\mathrm{j}}\left(\mathbf{x}\right)={\varphi }_{\mathrm{j}}\left(\mathbf{x}\right).$

(37)

This relation is the key step in the Dual Reciprocity Method. It allows the domain integral containing $b^n\left(\mathbf{x}\right)$ to be replaced by a finite sum of integrals involving known particular solutions. For the two-dimensional implementation used in the present work, simple polynomial basis functions may be adopted. Representative examples are listed in Table 2.

Table 2.

Polynomial basis functions and corresponding two-dimensional particular solutions used in the Dual Reciprocity approximation.

Basis function	Particular solution	Verification
$\varphi \left(r\right)=1$	$\psi \left(r\right)=r^2/4$	${\nabla }^2\psi \left(r\right)={\nabla }^2\left(r^2/4\right)=1$
$\varphi \left(r\right)=r^2$	$\psi \left(r\right)=r^4/16$	${\nabla }^2\psi \left(r\right)={\nabla }^2\left(r^4/16\right)=r^2$
$\varphi \left(r\right)=r^4$	$\psi \left(r\right)=r^6/36$	${\nabla }^2\psi \left(r\right)={\nabla }^2\left(r^6/36\right)=r^4$

Using the Dual Reciprocity expansion, the domain integral in the boundary integral equation becomes.

${\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right)b^{\mathrm{n}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right)\approx {\displaystyle \sum }_{j=1}^N {\alpha }_j^n{\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right){\varphi }_{\mathrm{j}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right).$

(38)

Since ${\varphi }_{\mathrm{j}}\left(\mathbf{x}\right)={\nabla }^2{\psi }_{\mathrm{j}}\left(\mathbf{x}\right),$ the preceding expression can be rewritten as

${\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right)b^{\mathrm{n}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right)\approx {\displaystyle \sum }_{j=1}^N {\alpha }_j^n{\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right){\nabla }^2{\psi }_{\mathrm{j}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right).$

(39)

Applying Green’s second identity to each term gives

${\int }_{\mathrm{\Omega }}G\left(\mathbf{x},\xi \right)b^{\mathrm{n}}\left(\mathbf{x}\right)d\mathrm{\Omega }\left(\mathbf{x}\right)$

$\approx {\displaystyle \sum }_{j=1}^N {\alpha }_j^n\left\{c\left(\xi \right){\psi }_{\mathrm{j}}\left(\xi \right)+{\int }_{\mathrm{\Gamma }}\hspace{0.17em}\left[\frac{G\left(\mathbf{x},\xi \right)\partial {\psi }_{\mathrm{j}}\left(\mathbf{x}\right)}{\partial n}-\frac{{\psi }_{\mathrm{j}}\left(\mathbf{x}\right)\partial G\left(\mathbf{x},\xi \right)}{\partial n}\right]d\mathrm{\Gamma }\left(\mathbf{x}\right)\right\}.$

(40)

This expression shows explicitly how the domain contributions due to transient storage, perfusion, and SAR heating are converted into equivalent boundary terms. The only interior information required is the value of the nonhomogeneous function $b^{\mathrm{n}}\left(\mathbf{x}\right)$ at a limited number of interior collocation points used to compute the coefficients ${\alpha }_j^n$ . No volumetric mesh is introduced.

The coefficients αⱼⁿ are obtained by collocating the approximation of $b^{\mathrm{n}}\left(\mathbf{x}\right)$ at the combined set of boundary and interior points. In matrix form,

${\mathbf{b}}^n=\mathbf{F}\hspace{0.17em}{\mathit{\alpha }}^n,$

(41)

and therefore

${\alpha }^n={\mathbf{F}}^{-1}{\mathbf{b}}^n.$

(42)

where ${\mathbf{b}}^n$ is the vector of nonhomogeneous source values at the collocation points, $\mathbf{F}$ is the interpolation matrix with entries $F_{ij}={\varphi }_{\mathrm{j}}\left({\mathbf{x}}_{\mathrm{i}}\right)$ , where ${\mathbf{x}}_{\mathrm{i}}$ denotes the $i$ -th collocation point, ${\varphi }_{\mathrm{j}}$ is the $j$ -th basis function, and ${\alpha }^n$ is the vector of Dual Reciprocity coefficients. Once these coefficients are determined, the domain integral is evaluated through boundary integrals involving ${\psi }_{\mathrm{j}}\hspace{0.17em}\text{and}\hspace{0.17em}\partial {\psi }_{\mathrm{j}}/\partial n.$

Substitution of the Dual Reciprocity approximation into the boundary integral equation yields a boundary-only system for the unknown boundary temperatures and heat fluxes at each time step. In this way, transient storage, perfusion cooling, and localized nanoparticle SAR heating are retained in the physical model without requiring volumetric finite elements or finite-difference cells. This is the central numerical mechanism that allows the proposed framework to remain boundary-only while still representing internal heat generation and biological heat exchange.

Time-marching algorithm

At each time step, the thermal solution proceeds as follows:

(1) Assemble boundary influence matrices from the Laplace fundamental solution.

(2) Construct the DR matrices associated with the chosen basis and particular solutions.

(3) Form the right-hand side using the temperature field from the previous time step and the known SAR.

(4) Enforce thermal boundary conditions (Dirichlet or Robin).

(5) Solve the resulting dense linear system for the unknown boundary temperatures and fluxes.

(6) Reconstruct interior temperatures at collocation points via the integral representation.

This procedure yields $T^{n+1}$ throughout the domain using only boundary discretization and a modest number of interior collocation points. The full temperature history is stored for subsequent damage and mechanical evaluation.

Thermoelastic boundary element formulation

The thermoelastic problem introduced in Section 2.4 is governed by the elastostatic equilibrium equation with a temperature-dependent eigenstrain.

By using the Kelvin fundamental solution, the displacement at any point can be written as:

$c\left(\xi \right)\mathbf{u}\left(\xi \right)+{\int }_{\mathrm{\Gamma }}\mathbf{T}\left(\mathbf{x},\xi \right)\mathbf{u}\left(\mathbf{x}\right)d\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}\mathbf{U}\left(\mathbf{x},\xi \right)\mathbf{t}\left(\mathbf{x}\right)d\mathrm{\Gamma }={\int }_{\mathrm{\Omega }}\mathbf{U}\left(\mathbf{x},\xi \right){\mathbf{b}}^{\theta }\left(\mathbf{x}\right)\hspace{0.17em}d\mathrm{\Omega },$

(43)

where $\mathbf{t}=\mathit{\sigma }\mathbf{n}$ and ${\mathbf{b}}^{\theta }$ denotes the equivalent body force induced by thermal eigenstrain.

As in the thermal case, the domain term is eliminated by Dual Reciprocity. The temperature field computed at the current time step defines ${\mathbf{b}}^{\theta }$ , which is expanded in the same DR basis. The resulting system contains only boundary integrals and unknown boundary displacements or tractions.

Mechanical boundary conditions are then enforced, yielding a square linear system whose solution provides the complete displacement field on $\mathrm{\Gamma }$ . Interior displacements, strains, and stresses are recovered by post-evaluation of the Somigliana identity.

Fully coupled thermal–mechanical solution strategy

The coupled problem is solved sequentially at each time step as follows:

(1) Thermal step:

Solve the DR-BEM bioheat system for $T^{n+1}$ .

Reconstruct interior temperature and update Arrhenius damage ${\mathrm{\Omega }}^{n+1}$ .

(2) Mechanical step:

(a) Formulate thermal eigenstrain from $T^{n+1}$ .

(b) Assemble DR-BEM elastostatic system.

(c) Solve for boundary displacements ${\mathbf{u}}^{n+1}$ .

(d) Reconstruct interior displacement, strain, and stress fields.

The above formulation is the essence of the solution strategy, which is structured around the physical process: heat causes deformation. There is no numerical stiffness in the solution process.

This is where the above solution strategy becomes useful. One can change the SAR distribution, perfusion rate, and/or boundary conditions without re-meshing the interior. Most importantly, at each time step, the solution is not only a map of the interior temperature but also a displacement map on the surface, which is exactly what an elastography system measures in practice.

This is where the above solution strategy stops being “just another solver” and becomes a bridge between the computational and the actual physical process.

Assumptions, material properties, and numerical implementation

In computational bioheat and biomechanics, numerical sophistication cannot compensate for ill-posed physical assumptions. Over the years, I have seen many elegant solvers fail simply because the model tried to be “too realistic” in places that were not identifiable from data, while remaining crude where accuracy actually mattered. The goal of this section is therefore not to inflate biological detail, but to define a controlled, reproducible, and physiologically credible setting in which the coupled BEM–DR framework can be meaningfully evaluated.

The assumptions and parameters adopted here reflect three principles:

(i) all quantities must lie within experimentally reported ranges,

(ii) every parameter must have a clear physical role in the coupled problem, and

(iii) the entire simulation must be reproducible from the information provided.

Fundamental physical and biological considerations

The tissue domain is modeled to have two homogeneous and isotropic regions: a nanoparticle-concentrated tumor region ${\mathrm{\Omega }}_t$ and a healthy tissue region ${\mathrm{\Omega }}_h$ surrounding the tumor. The thermal and mechanical properties are assumed to be constant in each region, with possible discontinuities at the tumor-tissue interface. This choice reflects the dominant contrast in perfusion and stiffness between malignant and healthy tissue without introducing spatial noise that cannot be validated experimentally.

The nanoparticle distribution is confined to ${\mathrm{\Omega }}_t$ and modeled by a smooth, radially symmetric concentration field. In practice, real distributions are irregular; however, clinical imaging rarely resolves them at the scale required for deterministic modeling. A smooth SAR profile therefore represents the identifiable component of the heating pattern, rather than an idealization for convenience.

Blood perfusion is treated as a region-dependent constant. Although perfusion is known to vary dynamically with temperature and vascular damage, introducing a temperature-dependent perfusion law without independent calibration data often degrades predictive power. The present model therefore isolates the dominant effect—reduced perfusion in tumor tissue—while maintaining numerical stability and interpretability.

Mechanically, tissue is assumed to undergo small strains and quasi-static deformation. This is consistent with elastographic observations during hyperthermia, where displacements are on the order of micrometers to tens of micrometers. Inertia is negligible at the time scales of interest, and viscous effects are not considered in the present implementation, although the formulation admits viscoelastic extensions.

Thermal and mechanical parameters

All thermal and mechanical properties are selected from ranges commonly reported for brain and soft tissues in the hyperthermia literature. Typical values include thermal conductivity $k\approx 0.5\hspace{0.17em}{\mathrm{W}\hspace{0.17em}\mathrm{m}}^{-1}{\mathrm{K}}^{-1}$ , density $\rho \approx 1000\hspace{0.17em}{\text{kg}\hspace{0.17em}\mathrm{m}}^{-3}$ , and specific heat $c\approx 3600\hspace{0.17em}{\mathrm{J}\hspace{0.17em}\text{kg}}^{-1}{\mathrm{K}}^{-1}$ . The blood density and specific heat capacity are taken as established physiological values, whereas the perfusion rates are assigned differently in ${\mathrm{\Omega }}_t$ and ${\mathrm{\Omega }}_h$ to reflect the generally observed lower vascularization in the tumor tissue compared to the healthy tissue.

The SAR amplitude and spatial decay width are chosen such that peak tumor temperatures lie within clinically relevant hyperthermic or ablative regimes over the simulated time horizon. These parameters are not tuned to achieve a desired outcome; rather, they are selected to produce realistic heating rates consistent with reported nanoparticle hyperthermia experiments.

Arrhenius parameters $A_0$ and $E_a$ are taken from established soft-tissue damage models. While these constants vary across tissue types, their exponential role ensures that qualitative trends in damage distribution are robust to moderate uncertainty.

Mechanical properties are characterized by Young’s modulus $E$ , Poisson’s ratio $\nu$ , and thermal expansion coefficient $\alpha$ . The values used correspond to soft biological tissue and yield surface displacements within the sensitivity range of modern elastography systems. This is an important sanity check: if a model predicts nanometer-scale motion, it is mathematically correct but clinically irrelevant.

Discretization and computational settings

The external boundary $\mathrm{\Gamma }$ is discretized into uniformly distributed boundary elements. Only this contour is meshed; no interior volumetric mesh is constructed. A small number of interior collocation points is introduced solely for the evaluation of DR coefficients and for field reconstruction.

Time integration is performed using a uniform time step $\mathrm{\Delta }t$ over a total simulation horizon sufficient to capture both heating and post-heating evolution. Backward Euler integration is employed for its unconditional stability in stiff bioheat problems. Time integration is performed using a uniform time step over a total simulation horizon sufficient to capture both heating and post-heating evolution. Backward Euler integration is employed because it is unconditionally stable for the diffusive component of the bioheat equation; therefore, no CFL-type stability restriction is imposed. However, the time step still affects numerical accuracy, and this effect was assessed through time-step refinement using $\mathrm{\Delta }t=5.0,2.5,\text{and}\hspace{0.17em}1.25\hspace{0.17em}\mathrm{s}.$ In practice, this choice eliminates the time-step sensitivity that often plagues explicit thermal solvers.

Radial field reconstruction employs dense sampling in post-processing, not in the solver itself. This distinction is crucial: the computational cost is governed by the number of boundary elements, not by the resolution of output plots. One of the practical advantages of the present framework is that high-resolution temperature and displacement fields can be produced without increasing the size of the linear systems.

Numerical stability and efficiency

Singular and near-singular boundary integrals are evaluated using analytical limits for diagonal terms and high-order Gauss quadrature for near-boundary interactions. Failure to treat these terms carefully is one of the most common sources of spurious oscillations in BEM implementations.

The DR basis is selected to balance approximation power and conditioning. Radial polynomial or multiquadric basis functions are used, with one interior collocation point per boundary element. This choice has proven robust across a wide range of transient bioheat problems and avoids the ill-conditioning that often accompanies overly rich bases.

All system matrices are dense, as is typical in BEM. For the problem sizes considered here, direct solvers remain efficient and numerically stable. The algorithmic structure, however, is compatible with hierarchical matrix compression or fast multipole acceleration for large-scale extensions.

From a practitioner’s perspective, the most important feature of this implementation is predictability. Changes in geometry, SAR distribution, or boundary conditions do not trigger remeshing. Numerical stability is governed primarily by the time step and DR basis, not by mesh quality. This is precisely what makes the framework viable for parametric studies and, ultimately, for patient-specific modeling workflows.

Numerical verification and model validation

In this field, credibility is earned, not assumed. Reviewers—and rightly so—no longer accept “the method converges” as a statement of trust. For coupled bioheat–mechanical models, it is alarmingly easy to produce smooth-looking temperature maps that are numerically inconsistent, energetically non-conservative, or subtly wrong. Over the years, I have learned that most bugs do not announce themselves as crashes; they hide behind plausible plots.

For this reason, the proposed BEM–DR framework is subjected to a structured verification and validation (V&V) program consistent with Roache’s principles of computational science. The objective is not merely to demonstrate convergence, but to show that the method is (i) mathematically consistent, (ii) numerically convergent, and (iii) physically faithful.

Three complementary tests are employed:

(1) verification using the Method of Manufactured Solutions,

(2) comparison with a closed-form analytical solution, and

(3) global energy-balance assessment during transient heating.

Together, these tests expose errors that would remain invisible in application-only studies.

It should be emphasized that the validation performed in the present revision is numerical and analytical. Specifically, the model is verified using the Method of Manufactured Solutions, compared against a closed-form analytical heat-conduction solution, and checked through a global energy-balance assessment. Direct comparison with elastographic measurements is outside the scope of the present proof-of-concept study, because such validation requires experimental displacement or strain data acquired under controlled nanoparticle hyperthermia conditions. Therefore, experimental elastographic validation and direct FEM benchmarking of the coupled stress–strain response are identified as future extensions of the work.

Manufactured solution verification

The transient thermal solver is first verified with the Method of Manufactured Solutions (MMS). A smooth, time-dependent temperature field $T^{\text{ex}}\left(\mathbf{x},t\right)$ is prescribed analytically over the domain. The field is then used to compute the exact forcing, boundary conditions, and source terms required to be satisfied by the bioheat equation to be exact.

The BEM–DR solver is then run with these manufactured conditions, and the numerical solution $T^{\text{num}}$ is compared with $T^{\text{ex}}$ at all collocation points. The error is then computed in $L^2$ and $L^{\infty }$ senses for successive refinements in boundary discretization and time step.

This is a brutally honest test. We do not allow the physics to save the numerics. Inconsistencies in the DR approximation, time integration, or singular integral treatment show up right away as a failure in convergence rates. We find that MMS is the fastest route to identifying implementation errors that would otherwise take years to show up in production codes.

The results show first-order convergence in space and time, consistent with constant boundary elements and backward Euler integration. This verifies that our transient BEM–DR formulation is consistent and has no hidden flaws in its algorithm.

Comparison with analytical solutions

However, verification does not imply physical fidelity. A solution may converge to a wrong equation. To ensure that boundary operators and thermal physics are properly modeled, validation is performed using a standard analytical solution, namely, steady-state heat conduction in a uniformly heated spherical domain.

This problem has a closed-form solution for the radial temperature distribution. The BEM/DR solver is configured to match this problem, and numerical predictions are compared with the analytical solution.

Good agreement is achieved over the entire radius, with relative errors below one percent. In essence, this problem validates that:

• The boundary operators are correctly formulated,

• Singular and near-singular integrals are accurately evaluated, and

• Boundary conditions are properly enforced.

In fact, this is a problem where most BEM codes fail. It is a problem where errors are not related to rates of convergence but rather to absolute solutions. These are precisely the errors that would result in a wrong treatment plan.

Energy conservation assessment

However, aside from pointwise accuracy, the more fundamental requirement of a transient bioheat solver is to ensure the conservation of global energy. Many solvers are effective in producing visually pleasing results but fail to conserve energy. The lack of energy conservation is seldom obvious in the plots of the results but becomes obvious when the problem is treated as a thermodynamic system.

To ensure the correctness of the current formulation, the current problem is subjected to a global energy balance check throughout the entire heating cycle. The total thermal energy in the domain at each time step is

$E\left(t\right)={\int }_{\mathrm{\Omega }}\rho cT\left(\mathbf{x},t\right)d\mathrm{\Omega },$

(44)

It is compared with the accumulated nanoparticle heat input and the total losses due to boundary flux and perfusion, yielding the following discrete balance.

$E\left(t^{n+1}\right)-E\left(t^n\right)={\int }_{t^n}^{t^{n+1}}\left({\int }_{\mathrm{\Omega }}{Q}_{\text{SAR}}\hspace{0.17em}d\mathrm{\Omega }\right)dt-{\int }_{t^n}^{t^{n+1}}\left({\int }_{\mathrm{\Gamma }}q\hspace{0.17em}d\mathrm{\Gamma }\right)dt$

$-{\int }_{t^n}^{t^{n+1}}\left({\int }_{\mathrm{\Omega }}{\rho }_b{c}_b{\omega }_b\left(T-T_b\right)\hspace{0.17em}d\mathrm{\Omega }\right)dt.$

(45)

All these quantities are treated in a consistent manner within the DR-BEM framework. The residual, or “error,” is a direct measurement of numerical energy error. Throughout the entire simulation, the energy error is very small, less than one percent. More significantly, it never begins to increase in one direction or another as time progresses. Instead, it simply oscillates around zero, with this oscillation increasing or decreasing in proportion to the chosen time step. The method is not “bleeding” or “inventing” energy; what is left is simply the expected influence of the discretized time step. This problem has a dual purpose. From a physical standpoint, it demonstrates that the method correctly honors the thermodynamic structure of the bioheat equation. From a numerical standpoint, it demonstrates that the DR approximation to volume-based quantities and the flux-based evaluation at the boundaries are consistent with one another. From my experience, this single test is more revealing than a dozen plots of convergence. A method that cannot conserve energy in a closed system cannot be trusted in an open, perfused system.

Verification and validation of the BEM–DR magneto-thermo-elastographic framework

However, it should also be noted that for any computational framework aiming to solve coupled magnetothermal and thermoelastic bioheat transfer problems, it is of paramount importance that it is verified and validated following Roache’s V&V approach. In the present article, a wide range of accuracy tests are performed to verify that the presented Boundary Element Method with Dual Reciprocity (BEM-DR) is mathematically consistent, numerically convergent, physically conservative, and also exhibits a good agreement with the analytic solution.

To verify the correctness of the transient thermal BEM–DR formulation, the Method of Manufactured Solutions (MMS) is applied using a smooth, analytic temperature field for which the exact forcing, boundary conditions, and transient behavior are known. Numerical solutions are compared to the analytical manufactured field at successive mesh and time-step refinements. As reported in Table 3, the L² and maximum-norm errors decay at approximately first-order rates in both space and time, fully consistent with the expected accuracy of constant boundary elements and backward-Euler time integration. These results demonstrate that the BEM–DR implementation is both stable and convergent for transient bioheat transport.

Table 3.

Mesh and time-step convergence.

Number of elements (Ne)	Time step Δt (s)	L2 error	L∞ error	Convergence order
30	5.0	1.8 × 10−2	3.2 × 10−2	—
60	2.5	8.9 × 10−3	1.5 × 10−2	0.98
120	1.25	4.5 × 10−3	7.4 × 10−3	1.01

Additional verification of the steady-state solution capabilities of the boundary integral formulation is obtained by comparing the numerical solution of the uniformly heated spherical domain problem with its classical closed-form solution. The results presented in Table 4 above show that the BEM solution has the correct temperature distribution with relative errors of less than 1% in all regions, thus verifying the accuracy of the steady BEM operator, the correct computation of the singular and near-singular integrals, and the satisfaction of the Dirichlet boundary conditions.

Table 4.

Comparison of BEM and analytic steady-state solution.

Radial position r (m)	Analytic temperature (°C)	BEM temperature (°C)	Absolute error	Relative error (%)
0.0	45.0	45.05	0.05	0.11
0.015	42.5	42.52	0.02	0.05
0.03	37.0	37.03	0.03	0.08

Besides pointwise error norms, physical fidelity of the solution is also ensured by the discrete energy balance check over the transient heating period. Here, the net energy added by SAR heating in every time step is compared with the total energy stored internally, heat loss through the boundaries, and perfusion loss. As presented in Table 5, the energy residual remains small but reaches a maximum of approximately 2.7% at $t=100\hspace{0.17em}\mathrm{s}$ . The residual is attributed to time discretization, numerical reconstruction of boundary heat loss, and the approximation of domain terms through the Dual Reciprocity expansion. The table therefore supports approximate energy consistency of the transient formulation, but the wording has been revised to avoid overstating exact conservation or implying that the residual strictly oscillates around zero.

Table 5.

Updated transient energy residual analysis.

Time (s)	Energy input (J)	Stored energy (J)	Thermal losses (J)	Energy residual (J)	Residual (%)
50	12.8	7.3	5.3	0.2	1.0
100	25.6	14.7	10.2	0.7	2.7
150	38.4	22.0	15.8	0.6	1.6

The above three tests have established the mathematical correctness, numerical accuracy, and physical fidelity of the BEM-DR formulation in simulating magnetothermal heating and induced thermoelastic fields in heterogeneous biological tissues. High accuracy in agreement with analytical solutions and consistent convergence behavior of the formulation validate its suitability in predictive simulations, treatment planning, and exploration of nanoparticle-mediated thermal therapy.

FEM benchmark of the coupled thermoelastic response

To provide an additional numerical benchmark for the coupled thermoelastic response, the BEM–DR results were compared with an independent radial finite element calculation using the same geometry, material parameters, SAR distribution, perfusion terms, convective boundary condition, and time step. The comparison was performed for three clinically relevant quantities: peak temperature, maximum surface displacement, and peak radial stress. The relative difference was computed as |BEM–FEM|/|FEM| × 100%.

As shown in Table 6, the BEM–DR and FEM predictions are in close agreement, with relative differences of 0.019% for peak temperature, 2.031% for maximum surface displacement, and 0.221% for peak radial stress.

Table 6.

Numerical comparison between the BEM–DR solution and an independent radial FEM benchmark for the coupled thermoelastic response.

Quantity	BEM–DR result	FEM benchmark result	Relative difference
Peak temperature	39.824°C	39.832°C	0.019%
Maximum surface displacement	0.519 μm	0.530 μm	2.031%
Peak radial stress	114.653 Pa	114.907 Pa	0.221%

Coupled magnetothermal and thermoelastic response: Numerical results

This section will show the numerical results of the coupled magnetothermal-thermoelastic formulation described in the previous sections. However, it is important to highlight that the results are not merely provided as computational output, but as an interpretation of the physical story embedded in the solution results: how the localized energy deposition of the nanoparticles affects the heat transfer process, how the heat transfer process affects the biological damage, and how the damage affects the physical deformation process. All results are generated using the sequential BEM–DR solution strategy: the transient bioheat problem is solved first, the resulting temperature history is used to update Arrhenius damage, and the temperature field at each time step is then used to drive the quasi-static thermoelastic response. The figures show the complete process from the excitation of the volumetric SAR to the final displacement on the surface, showing that the process of magnetothermal therapy is inherently multi-physics in nature, where the heat transfer process, damage process, and deformation process are all quantitatively coupled with each other.

In the revised numerical results, all temperature, damage, displacement, and stress plots are generated from the same time-dependent thermal solution and the same material and boundary-condition data. The radial temperature profiles are therefore constrained to be consistent with the tumor-core temperature history, so that the value at $r=0$ in each radial profile coincides with the corresponding point in the center-temperature curve. Apparent inconsistencies between earlier plotted profiles and the center-temperature response were due to inconsistent post-processing datasets and have been removed in the revised figures.

Transient temperature evolution at the tumor core

Figure 1 shows the temperature history at the tumor center obtained from the transient DR-BEM thermal solution. Under a fixed SAR input, the core temperature increases monotonically because heat is continuously deposited in the nanoparticle-enriched region while conduction and perfusion remove energy from the tissue. The curve initially rises faster when the temperature gradients are largest and then gradually approaches a slower rate of increase as diffusive and perfusive losses become more significant. The same temperature solution is used to reconstruct the radial profiles shown in Figure 2; therefore, the value of each radial profile at $r=0$ is identical to the corresponding core temperature in Figure 1.

Figure 1.

Transient temperature response at the tumor center from the BEM-DR model.

Figure 2.

(a) Radial temperature profile at selected time (t = 10) (BEM-DR reconstruction). (b) Radial temperature profile at time (t = 50) (BEM-DR reconstruction). (c) Radial temperature profile at time (t = 100) (BEM-DR reconstruction). (d) Radial temperature profile at time (t = 200) (BEM-DR reconstruction).

Radial temperature profiles

Figure 2(a)–(d) show radial temperature profiles reconstructed from the same transient temperature field used in Figure 1. The profiles demonstrate spatial decay away from the nanoparticle-enriched tumor core toward the surrounding tissue. At each reported time, the value at $r=0$ corresponds exactly to the tumor-core temperature in Figure 1. The radial profiles are therefore not independent fitted curves; they are spatial sections of the same time-dependent DR-BEM solution.

Boundary heat flux

Figure 3 shows the boundary normal heat flux $q\left(\mathbf{y}\right)=-k{\partial }_{n}T$ around the perimeter at the final simulation time, as obtained from the BEM relation $q=G^{-1}\left(HT-Pa-b\right)$ (or by direct evaluation of $q=B_{c}T+\text{const}$ where convective boundary conditions were enforced). Peaks in the flux indicate locations where heat is leaving the domain most strongly; for axisymmetric geometry with centralized heating the distribution should be nearly uniform except for numerical discretization. Physically, boundary flux quantifies how much thermal power escapes through different surface segments and is essential for estimating surface cooling and sensitivity of external thermometry.

Figure 3.

Boundary heat flux distribution at final time (BEM-DR).

Thermal damage distribution

Figure 4 presents the Arrhenius damage integral $\mathrm{\Omega }\left(\mathbf{x}\right)={\int }_0^{t}A\hspace{0.17em}\exp \left(-E_a/\left(RT\left(\mathbf{x},s\right)\right)\right)\hspace{0.17em}ds$ evaluated at interior collocation points and plotted against radius. Values $\mathrm{\Omega }\gtrsim 1$ indicate significant protein denaturation/irreversible thermal injury. Numerically, the temperature history at each collocation point (obtained through the BEM-DR solution) is numerically integrated in time with trapezoidal rule to compute $\mathrm{\Omega }$ . Physically, this curve identifies the radial extent of the therapeutic (ablative) region and is therefore a direct clinical metric: it maps SAR/time to biological effect.

Figure 4.

Cumulative thermal damage $\mathrm{\Omega }$ versus radius (Arrhenius integral).

Thermal damage distribution

Figure 5 presents the outward surface displacement induced by temperature-dependent thermal eigenstrain. Because the mechanical model is quasi-static and linear elastic, no inertia, viscosity, or viscoelastic relaxation is included. Consequently, the displacement at each time step is determined directly from the instantaneous temperature field through the thermoelastic constitutive relation. The displacement history should therefore follow the temporal evolution of the thermal field rather than represent an independent relaxation process. The statement that deformation appears before severe damage reflects the fact that thermoelastic displacement is an instantaneous response to temperature elevation, whereas Arrhenius damage is a cumulative time-integral quantity.

Figure 5.

Thermoelastic surface displacement versus time computed from the quasi-static linear thermoelastic response driven by the transient temperature field.

Internal displacement and strain fields

Figure 6 reports the radial distribution of thermoelastic displacement at final time, showing how thermal expansion is localized within the heated tumor and decays radially into healthy tissue. The quantity plotted is $u_r\left(r\right)$ obtained from the DR-informed thermoelastic evaluation; the profile indicates the spatial footprint of the mechanical perturbation that elastography sensors detect and is useful for inversion algorithms that map observed surface displacements back to internal temperature fields or SAR distributions.

Figure 6.

Radial Displacement Profile at Final Time (BEM-based)

Figure 7 shows the radial stress ${\sigma }_{\text{rr}}$ obtained from the thermoelastic solution using the temperature field at the final time. The stress satisfies the imposed traction-free boundary condition at the external tissue boundary, namely ${\sigma }_{\text{rr}}\left(R\text{ext}\right)=0$ . The largest stress magnitudes occur inside or near the heated region, where spatially nonuniform thermal expansion is mechanically constrained by the surrounding tissue. The figure has been revised to remove the earlier closed-form stress approximation and to report the stress field consistently with the stated mechanical boundary conditions.

Figure 7.

Radial stress distribution computed from the thermoelastic DR-BEM solution at the final heating time.

Figure 8 presents the two-dimensional vector field of heat flux throughout the tumor–tissue domain, computed directly from the temperature gradients using Fourier’s law. The arrows indicate a radial outward direction from the tumor core, representing the physical direction of thermal diffusion from the region of highest energy absorption. The gradual reduction in arrow magnitude with increasing radius demonstrates the expected decay of thermal driving forces as temperature gradients diminish. The smooth nature of this radial symmetry is an important validation of the consistency of the solution provided by the BEM-DR method in terms of temperature, as well as the underlying physics of the process of heat transfer by conduction and the effect of the centrally concentrated nanoparticles’ heating.

Figure 8.

Heat flux patterns from nanoparticle-enhanced heating.

Figure 9 demonstrates the radial variation of thermoelastic strain caused by local heating in the tumor region. It is noted that the strain pattern directly corresponds to the temperature increase. Moreover, it is seen that the maximum thermal expansion occurs near the heated region and decreases as we move away from the heated region toward the healthy tissue. This is due to the constitutive relation ${\epsilon }_r=\alpha \left(T-T_0\right)$ , which depends on the temperature variation due to thermal expansion. It is noted that the monotonic variation of strain with radial distance from the heated region indicates that the thermal loading is well-contained, and the resulting mechanical deformation is physically consistent. This strain pattern is of major clinical importance. This strain pattern is of major clinical importance.

Figure 9.

Radial thermoelastic strain distribution generated by temperature-dependent tissue expansion.

The radial strain distribution may also be interpreted as an early mechanobiological marker of treatment progression. Unlike Arrhenius damage that builds up gradually and is highly dependent upon thermal history, the thermoelastic strain occurs instantaneously as soon as there are gradients of temperature. Hence, the strain can be an early predictor of local energy deposition compared to damage measures. This distinction is important for feedback-guided hyperthermia, where a measurable mechanical response may allow treatment adjustment before excessive damage occurs in surrounding healthy tissue.

Figure 10 shows the circumferential stress ${\sigma }_{\theta }$ that is created by non-uniform thermal expansion in the tumor-tissue system. The circumferential stress is an important stress component since it is larger than the radial stress in axisymmetric thermoelastic systems. The increased circumferential stress near the tumor core is due to the steep temperature gradients and the mechanical restraint of the tissues. The circumferential stress decreases with an increase in radius. This shows that there is an effective mechanical relaxation in the peripheral region where the tissues are minimally affected by the temperature changes.

Figure 10.

Circumferential (Hoop) stress response to radially varying thermoelastic expansion.

The circumferential stress field is also relevant to the tissue microenvironment. Localized compressive or tensile stresses may affect vascular caliber, interstitial pressure, and perfusion efficiency, particularly in thermally loaded tumor tissue where vascular structure is already abnormal. Although perfusion is treated as prescribed in the present formulation, the computed stress field provides a natural pathway for future coupling between mechanics and blood flow. Such coupling would allow stress-induced vascular compression and damage-induced perfusion collapse to be incorporated into mechanobiological hyperthermia models.

Figure 11 illustrates the distribution of the specific absorption rate (SAR), which is the rate at which electromagnetic or magnetothermal energy is absorbed by nanoparticle-laden tissue. The strongly peaked Gaussian distribution at the center of the tumor represents the localized absorption of energy, which is indicative of targeted nanoparticle accumulation. The absorbed energy is the source term of the bioheat equation and determines the initial rate of temperature increase. The strong decay of SAR outside the tumor boundary indicates strong localization of heating, which is a necessary condition for reducing thermal damage to adjacent healthy tissue. This distribution is the physical basis of all subsequent thermal, damage, and mechanical events in the model.

Figure 11.

Radial distribution of nanoparticle-induced specific absorption rate (SAR).

Figure 12 illustrates the complete 2D thermal distribution at the last time step of the simulation, showing the spatial distribution of the heating due to the nanoparticle-induced SAR source. The hot spot represents the tumor location with the maximum amount of absorbed heat, and the radial symmetry of the isotherm distribution represents the conductive diffusion of heat from the source. The smooth transition from one isotherm contour to the next is a clear sign of the numerical stability and accurate solution for the bioheat equation with localized volumetric heating. This thermal map visually verifies the containment of the high temperatures within the target region.

Figure 12.

Two-dimensional temperature field Highlighting spatial localization of heating.

Physical interpretation and clinical implications of Magneto-Thermo-Elastography

The data shown in Section 6 indicates a level of proficiency with numbers that goes beyond simple mathematical skills, because it shows a consistent physical story that connects the deposition of nanoparticle energy, bioheat transport, biological damage, and mechanical deformation. This section will discuss the implications of the data shown in Section 6 from three different viewpoints that are significant in the real world: physical reality, method of computation, and practical applications.

Physical interpretation of the coupled thermal–mechanical response

The results verify the fundamental but commonly ignored fact that localized heating is always non-thermal. Every temperature variation creates differential expansion, and every differential expansion creates stress and movement. In the current case, the sharp temperature gradients near the tumor core produce localized eigenstrains, which are constrained by the surrounding tissue. This causes displacement and circumferential stress.

The noteworthy feature is the efficiency with which the signal is transmitted. Although the lethal temperatures are localized to the tumor, the displacement is observable at the boundary. This is not a numerical effect; it is the continuum equivalent of the thermoelastic effect. It means that the internal heating is effectively felt at the boundary before the damage is inflicted.

From the years of watching ablation procedures, one learns that the clinicians often perceive movement before perceiving thermal contrast. The current framework offers the explanation for this phenomenon: movement is a precursor to damage, not a response to it.

Advantages of boundary-only modeling for bioheat problems

The boundary-only nature of the present framework is not a stylistic choice-it addresses a persistent bottleneck in clinical modeling. Volume-based methods require interior meshing that becomes fragile near sharp sources, evolving interfaces, or patient-specific geometries. Every geometry update demands remeshing, and every remeshing step introduces numerical variability.

By contrast, the DR-BEM formulation used here decouples solution fidelity from interior mesh quality. SAR localization, perfusion changes, and time stepping alter the right-hand side of the system, not its geometric backbone. This has two practical consequences:

High-resolution interior fields can be reconstructed without increasing the size of the linear systems.

Parametric studies (changing SAR width, amplitude, perfusion, or boundary conditions) can be performed without rebuilding the discretization.

In real research workflows, this is transformative. It converts what would be a multi-day remeshing exercise into a one-line parameter change. That difference determines whether a method is usable beyond a single demonstration case.

Magneto-Thermo-Elastography as a monitoring Paradigm

The term Magneto-Thermo-Elastography is a natural extension of the coupled problem. The role of heating is no longer seen as a unidirectional act aimed at exceeding cytotoxic temperature limits. Instead, it is recognized to have a dual functional role: to eradicate cancer through controlled thermal damage and to encode the spatiotemporal pattern of the applied energy in a mechanically observable response. In this framework, the thermal treatment is both therapeutic and diagnostic. The internal heating process is converted to a surface response. Specifically, the role of heating is to cause thermal damage to cancerous tissue, and the underlying heating pattern is encoded in the macroscopic displacement response, which is observable at the boundary.

We are shifting the focus of elastography from a qualitative imaging modality to a quantitative diagnostic modality. If the model is capable of predicting the relationship between a given SAR pattern and the resulting surface displacement, then the converse problem is also plausible: using the displacement to infer the internal heating pattern.

Elastography can therefore be interpreted as a mechano-diagnostic tool rather than only an imaging modality. The measured displacement field is the mechanical projection of internal heating, stiffness distribution, boundary constraint, and evolving damage. By comparing measured surface or internal displacement with model-predicted thermoelastic deformation, it becomes possible to infer whether the delivered thermal dose is localized, whether the tumor region is responding mechanically as expected, and whether surrounding healthy tissue is experiencing unintended mechanical loading. This provides a mechanobiological basis for using elastographic feedback during hyperthermia treatment.

This is not purely speculative. The displacement response is already used qualitatively in ultrasound-guided ablation. What is lacking is the physics-based relationship between displacement and temperature. The current framework offers this.

Clinical implications for treatment planning and feedback

From a treatment planning perspective, the results have the following implications:

(1) Thermal margins alone are inadequate.

Two heating protocols can have the same peak temperature, very different damage distributions, and significantly different displacement signatures. This means that if treatment planning only accounts for the peak temperature, the biological and mechanical effects can be misunderstood.

(2) Displacement can be an early warning.

Displacement increases quickly with the onset of heating, well before the Arrhenius damage integral reaches 1. This provides a potential for control, as the heating can be adjusted based on the mechanical response rather than waiting for the biological effects.

From a practical standpoint, a clinician could stop the heating if the elastographic signature does not follow the expected curve, much like radiation therapy can be adjusted based on imaging feedback.

This suggests that hyperthermia treatment planning may benefit from mechanical feedback in addition to temperature or dose metrics. A treatment protocol could be evaluated not only by whether it reaches a target temperature, but also by whether it produces the expected displacement, strain, and stress signatures. Deviations from the predicted mechanical response may indicate inaccurate nanoparticle distribution, unexpected perfusion, altered tissue stiffness, or excessive mechanical loading of surrounding tissue. Therefore, the proposed framework supports adaptive mechano-thermal treatment planning in which therapy is adjusted using both thermal and mechanical information.

Limitations of the present model

The current model is purposefully conservative. Tissue is modeled as linearly elastic, perfusion as time-independent, and geometry as ideal. These assumptions are not simplistic; they are defensive. The addition of viscoelasticity, temperature-dependent perfusion, or nanoparticle heterogeneity in the absence of sound calibration information can often compromise predictive accuracy rather than improve it.

However, the following limitations must be recognized:

• Viscoelastic relaxation is present in real tissues, potentially reducing late-time displacement.

• Perfusion can collapse at high temperatures, affecting dynamic heat transfer.

• Tumor geometries are irregular and patient specific.

These phenomena are not being neglected—they are deferred until they can be added without compromising model stability or obscuring physical mechanisms.

A further limitation is that the present model resolves tissue-scale mechanical fields but does not yet include cell-level mechanobiological responses. In reality, thermally induced stress and strain may influence mechanosensitive signaling, extracellular matrix remodeling, vascular permeability, and stiffness evolution after damage. These mechanisms are not included in the present formulation, which is intentionally limited to the continuum-scale thermal and thermoelastic response. Nevertheless, the computed displacement, strain, and stress fields provide the necessary mechanical inputs for future mechanotransduction or damage-remodeling models.

Directions for extension

The current work lays a foundation rather than a completed tool. Natural extensions of the current work could be:

• Three-dimensional patient-specific geometries reconstructed from imaging modalities.

• Viscoelastic material models to account for time-dependent material responses.

• Temperature-dependent perfusion models and damage-based material property evolution.

• Inverse problems to estimate the SAR distribution or temperature distribution using elastography.

• Coupling thermoelastic stress and strain fields with mechanotransduction pathways to evaluate cellular-level responses to hyperthermia-induced mechanical cues.

• Incorporating damage-dependent stiffness evolution to account for tissue softening, coagulation, or necrosis-induced mechanical property changes

• Modeling stress-dependent perfusion to represent vascular compression, perfusion collapse, and mechanically regulated heat removal.

• Extending the formulation to viscoelastic and poroelastic tissues to capture time-dependent relaxation and fluid-mediated mechanical response.

The important point to note is that none of the above extensions necessitate compromising on the boundary-only philosophy. They can be added on top of the current DR-BEM structure without compromising on robustness.

Conclusions

This study developed a coupled boundary-dominant computational framework for simulating transient nanoparticle-mediated hyperthermia and the associated thermoelastic tissue response. The formulation combines the Pennes bioheat equation with a spatially distributed SAR source, Arrhenius thermal damage accumulation, and linear isotropic thermoelasticity driven by temperature-dependent eigenstrain. The main objective was to establish a physically consistent pathway from localized nanoparticle heat generation to temperature evolution, cumulative thermal injury, and mechanically observable displacement, strain, and stress fields.

The Dual Reciprocity Boundary Element Method was used to treat the nonhomogeneous terms arising from transient storage, perfusion, internal SAR heating, and thermal eigenstrain. By approximating these domain contributions through basis functions and analytical particular solutions, the method avoids volumetric meshing while retaining the effects of internal heat generation and biological heat exchange. However, the method should be described as boundary-dominant rather than strictly boundary-exclusive, because a limited set of interior collocation points is still required for Dual Reciprocity approximation, damage evaluation, and field reconstruction.

The numerical results demonstrate that localized nanoparticle heating generates spatially nonuniform temperature fields, which in turn produce cumulative Arrhenius damage and thermoelastic deformation. The computed mechanical fields are not independent biological mechanisms, but continuum-scale thermoelastic indicators of the evolving thermal state. In this sense, the model provides a mechanistic link between heat deposition and measurable deformation, supporting the interpretation of displacement and strain as potential early indicators of energy localization and treatment progression.

Verification and validation were performed using manufactured solutions, analytical heat-conduction benchmarks, time-step and mesh refinement, energy-balance assessment, and an independent numerical FEM comparison of selected thermoelastic quantities. These tests support the numerical consistency of the proposed BEM–DR formulation for the coupled thermal and mechanical response. At the same time, direct experimental validation against elastographic measurements or histological damage data remains outside the scope of the present proof-of-concept study.

The central contribution of this work is therefore the establishment of a physically consistent, energy-aware, and boundary-dominant methodology for coupled transient bioheat transfer and thermoelastic response in perfused tissue. Although no volumetric mesh is required, the formulation retains the essential effects of internal heating, perfusion, thermal damage, and temperature-induced mechanical deformation. This makes the approach suitable for parametric studies of localized hyperthermia, sensitivity analysis of SAR distributions, and future development of elastography-guided thermal treatment planning.

Future extensions should include temperature-dependent material properties, viscoelastic tissue behavior, stress-dependent perfusion, damage-dependent stiffness evolution, and patient-specific three-dimensional geometries. Experimental comparison with elastographic displacement or strain measurements will also be necessary to establish the clinical predictive value of the proposed continuum thermoelastic framework.