Introduction

High-Intensity Focused Ultrasound (HIFU) ablation is a non-invasive tumor treatment technology. The principle of HIFU is to focus ultrasound waves generated by the transducer onto a target tissue volume, inducing a rapid temperature rise. This localized heating causes coagulative necrosis of the tumor while sparing the surrounding and overlying healthy tissues [ 1 , 2 ]. HIFU thermal ablation has many advantages, such as non-invasive, non-contact, non-ionizing, and low-cost, and has been successfully used in the clinical treatment of solid malignant tumors, such as prostate cancer, liver cancer, kidney cancer, breast cancer, and pancreatic cancer [ 3 - 5 ].

The main idea of HIFU technology is to increase the temperature by focusing ultrasound power on the desired area to produce tissue coagulation or necrosis. Therefore, precise control of ultrasound distribution is significance. In previous studies, numerical simulations are often used to predict ultrasound propagation and thermal lesions caused by focused ultrasound [ 6-8 ].

Chen et al. [ 9 ] verified the nerve ablation effect of the HIFU transducer near the interface of bone and soft tissue through numerical simulation. The results showed that the use of transducers with incoherent modes can effectively optimize nerve ablation and reduce lesions to surrounding tissues. Huang et al. [ 10 ] studied and designed a dual-frequency transducer and compared the effects of single-frequency and dual-frequency ultrasound in thermal ablation. The results showed that dual-frequency switching can produce a larger ablation area and provide better control and flexibility. Zou et al. [ 11 ] proposed a coupled acoustic-thermal model to study the formation of thermal lesions and its optimal control in HIFU scanning treatment. The study analyzed the influence of adjacent treatment points on temperature rise and thermal lesion, and proposed a control scheme for optimizing heating time and sound intensity to improve treatment efficiency and uniformity. Kim et al. [ 12 ] proposed a numerical simulation tool coupling the Finite-Difference Time-Domain (FDTD) and Lattice Boltzmann (LBM) methods (FDTD-LBM) to predict the temperature distribution and lesion area of biological tissues in HIFU treatment. The study showed that this method has high accuracy in analyzing the thermal field caused by HIFU and can effectively predict the range of tissue lesions. Fura et al. [ 13 ] developed a numerical tool to predict the location and range of necrotic lesions in tissues after HIFU treatment. By comparing the numerical simulation results with the experimental data, the effectiveness of the model was verified, and the effects of different HIFU parameters were evaluated on the necrotic area. Fura et al. [ 14 ] explored the possibility of shortening the treatment time by optimizing HIFU treatment parameters, such as sound intensity and exposure time. Gharloghi et al. [ 15 ] optimized HIFU parameters for liver tumor ablation using MATLAB R2018b, focusing on frequency and power. The study found that 2 MHz and 100 W provided optimal thermal dose, reducing treatment time and damage to healthy tissues.

In previous studies [ 9-15 ], numerical simulations have been widely used to predict ultrasonic propagation and the formation of thermal lesions. Research has not fully explored how harmonic components, generated by acoustic wave propagation in nonlinear media, affect thermal lesions. While foundational studies, including our prior work [ 16 ], have demonstrated the value of harmonic monitoring, particularly the second harmonic ratio (P₂/P₁), in predicting thermal lesions under temperature-dependent parameters, they primarily captured the initial stage of nonlinear distortion. This study significantly advances the field by introducing a higher-order harmonic ratio, (P₃+P₄)/P₂, specifically designed to probe the efficiency of the subsequent nonlinear energy cascade into higher harmonics, which are more strongly absorbed and thus more directly responsible for localized heating. Furthermore, unlike the Spherical Beam Equation (SBE) model used in [ 16 ] at 1 MHz, we implement the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation optimized for deep ablation at 1.4 MHz, enabling a more accurate simulation of nonlinear propagation in liver tissue. Crucially, we established a dynamic, bidirectional acoustic-thermal coupling model that provides time-resolved evidence of how temperature-induced changes in attenuation (e.g., the reversal near 50 °C) directly govern this higher-order harmonic ratio, a mechanistic insight not previously reported. Finally, we rigorously compared thermal dose predictions using both the empirical Sapareto-Dewey and the mechanistic Arrhenius models, validating the applicability of our harmonic ratio biomarker across different lesion assessment frameworks. These combined methodological and analytical advancements provide a novel, more sensitive, and physically-grounded tool for real-time monitoring and control of HIFU therapy. This combination of a higher-order harmonic ratio biomarker, KZK-based modeling for deep tissue, dynamic acoustic-thermal coupling, and dual thermal dose validation provides a novel and more comprehensive tool for monitoring and controlling HIFU therapy, advancing beyond previous approaches that relied primarily on P₂/P₁ [ 16 ].

In particular, changes in the harmonic components generated during nonlinear propagation are critical. To quantify this, we defined a specific harmonic ratio:

$\begin{array}{cc}{R}^{*}=(P_3+P_4)/P_2& \mathrm{(1)}\end{array}$

where P_n is the pressure amplitude of the n-th harmonic. The selection of this formulation, (P₃+P₄)/P₂, is deliberately chosen over more common ratios (e.g., P₂/P₁) to create a more sensitive biomarker for the physical mechanisms driving thermal lesion formation. The rationale is twofold. First, the numerator, P₃+P₄, represents the cumulative pressure of higher harmonics that, due to their higher frequencies, are absorbed much more rapidly by tissue than the fundamental or second harmonic. This enhanced, localized energy deposition is a primary driver of rapid temperature rise and lesion creation. The second rationale, which is even more critical, is normalization by the second harmonic pressure, P₂. In the nonlinear energy cascade, the second harmonic acts as a powerful secondary source that generates the third and higher-order harmonics through wave-wave interactions. Through normalization by P₂, the ratio becomes largely independent of the absolute level of the second harmonic generated from the fundamental wave. Rather, it directly measures the conversion efficiency of P₂ energy into higher harmonics (P₃ and P₄). This makes the metric highly sensitive to temperature-induced changes in the tissue’s acoustic properties (e.g., attenuation and the nonlinearity parameter B/A), which directly govern the efficiency of this cascading process. Therefore, this ratio is not only a measure of overall nonlinearity, but a specific probe into the dynamics of the most thermally significant part of the acoustic energy spectrum. This formulation is supported by nonlinear acoustic theory, wherein the second harmonic serves as a source for higher harmonics via wave-wave interactions.

This study aimed to investigate the influence of harmonic ratio on HIFU-induced thermal lesion formation in biological tissue. By developing a coupled acoustic-thermal model, we simulate the propagation of HIFU waves and the resulting temperature elevation in tissue, incorporating nonlinear acoustic effects and dynamic tissue properties. The results of this study will offer significant contributions to enhancing HIFU therapy through the regulation of harmonic ratio, leading to improved treatment effectiveness and patient safety. Critically, nonlinear wave propagation generates harmonics that intensify heating, while temperature-dependent tissue properties (e.g., attenuation, sound speed) concurrently alter harmonic generation efficiency, a bidirectional acoustic-thermal coupling central to this study.

In summary, this study advances beyond prior work by introducing a higher-order harmonic ratio, employing a KZK-based model for deep tissue, establishing dynamic bidirectional coupling and validating with dual thermal dose criteria.

Material and Methods

In the numerical simulation of temperature rise and thermal tissue lesion, the widely used nonlinear KZK equation is used to describe the nonlinear ultrasonic propagation in tissues, which is expressed as [ 17 ]:

$\begin{array}{cc}\frac{{\partial }^{2}p}{\mathrm{\partial z}\partial \mathrm{t\prime }}=\frac{c_0}{2}{\Delta }_{\perp }p+\frac{b}{2{\rho }_0{c}_0^3}\frac{{\partial }^{3}p}{\partial {\mathrm{t\prime }}^2}+\frac{\beta }{2{\rho }_0{c}_0^3}\frac{{\partial }^2{p}^2}{\partial {\mathrm{t\prime }}^2}& \mathrm{(2)}\end{array}$

In Equation (2), the three terms on the right-hand side correspond to diffraction, attenuation, and nonlinear effects, respectively. The symbols p, ρ₀, c₀, b, and β represent the sound pressure, propagation medium density, sound wave velocity, dissipation parameter, and nonlinear coefficient, respectively. The Laplace operator of the rectangular coordinate system is ${\Delta }_{\perp }={\partial }^2/\partial x^2+{\partial }^2/\partial y^2,\mathrm{t\prime }=t-z/c_0$, represents the lag time, and z represents the coordinate along the acoustic axis.

Let $Z=z/F,R=r/a,\tau =\mathrm{\omega t},P=p/P_0$ and the normalized KZK equation is [ 17 ]:

$\begin{array}{cc}\frac{{\partial }^{2}P}{\mathrm{\partial \tau }\mathrm{\partial Z}}=\frac{1}{\mathrm{4G}}{\Delta }_{\perp }P+A\frac{{\partial }^{3}P}{\partial {\tau }^3}+\frac{N}{2}\frac{{\partial }^2{P}^2}{\partial {\tau }^2}& \mathrm{(3)}\end{array}$

where ω, P₀, a and F represent the angular frequency, the acoustic pressure on the transducer surface, the radius of transducer, and the focal length, respectively. $G\omega a^2/\left(2c_{0}F\right)$, $A=F\alpha {\omega }^2/\left(2c_0^3\right)$, where α is the attenuation coefficient at the fundamental frequency f₀. $N=F/l_d$, where $l_d={\rho }_0{c}_0^3/\left(\mathrm{\beta \omega }{P}_0\right)$ is the shock wave formation distance, describing the nonlinear effect, and ρ₀ represents the density. This equation can be combined with measurement calibration techniques to obtain boundary conditions and use the finite difference algorithm to calculate the solution in the frequency domain. In the frequency domain, the acoustic pressure is expressed in the form of a Fourier series expansion as follows [ 18 ]:

$\begin{array}{cc}p\left(\mathrm{z,\; r,\; t}\right)={\sum }_{\mathrm{-\infty }}^{\mathrm{+\infty }}{C}_n\left(\mathrm{z,\; r}\right)\exp \left(\mathrm{jn\omega \tau }\right)& \mathrm{(4)}\end{array}$

where $C_n\left(\mathrm{z,\; r}\right)$ represents the complex amplitude of the n-th harmonic component, and ω represents the fundamental frequency of the HIFU pulse.

Substituting equation (4) into equation (2), the amplitude of each harmonic can be obtained.

When ultrasound propagates in biological tissues, part of the sound energy is absorbed by the tissues and converted into heat energy. The temperature rise of the tissues is calculated using the Pennes bioheat equation [ 19 , 20 ].

$\begin{array}{cc}\frac{\mathrm{\partial T}}{\mathrm{\partial t}}=\frac{k_t}{{\rho }_0{C}_t}{\nabla }^{2}T+\frac{Q_v}{{\rho }_0{C}_t}& \mathrm{(5)}\end{array}$

Where C_t and k_t are tissue specific heat and thermal conductivity, respectively, T is tissue temperature, and Q_v is the amount of heat absorbed by the tissue per unit volume per unit time. The Pennes equation as implemented here does not include the blood perfusion cooling term, as this study focuses on establishing the fundamental relationship between harmonic ratio and thermal lesion formation in a controlled, ex vivo context. This simplification allows for the isolation of the acoustic-thermal-harmonic coupling mechanism from the effects of highly variable blood flow. This simplification is justified for ex vivo tissue but will be extended to in vivo models in future. Note that the blood perfusion term is omitted in this study to isolate the acoustic-thermal-harmonic coupling mechanism, consistent with an ex vivo experimental setup. Blood flow effects will be addressed in future tumor-specific models. Q_v can be calculated using the following formula [ 21 ]:

$\begin{array}{cc}{Q}_v=\frac{1}{{\rho }_0{c}_0}{\sum }_{\mathrm{n=1}}^N\mathrm{4\alpha }{n}^{\mu }{\left|C_n\right|}^2& \mathrm{(6)}\end{array}$

Here, N is the calculated harmonic number, μ is approximately 1 in gels and biological tissues.

According to the thermal dose formulation proposed by Sapareto and Dewey [ 22 , 23 ], the calculation of tissue lesion thermal dose is expressed as:

$\begin{array}{cc}\mathrm{CEM}\mathrm{43}={\int }_{t=t_0}^{t=t_1}{R}_T^{\mathrm{43}-T}\mathrm{dt}& \mathrm{(7)}\end{array}$

where CEM43 is the thermal dose received by the tissue. If T≥43 ℃, then R_T=0.5, otherwise R_T=0.25. This empirical formula expresses the thermal dose in Cumulative Equivalent Minutes at 43 °C (CEM43), which standardizes the biological impact of any time-temperature history. Experiments have shown that thermal necrosis will occur in tissues after being exposed to 43 ℃ for more than 240 minutes. This formula is also applicable to the judgment of thermal lesions in the study of HIFU treatment.

The thermal dose calculated using this empirical formula is expressed in Cumulative Equivalent Minutes (CEM), which quantifies the cumulative biological effect of a given thermal exposure by converting it into the equivalent number of minutes at the standard reference temperature of 43 ℃.

In addition, we implemented an Arrhenius injury integral to provide a mechanistic thermal lesion metric suitable for rapid heating [ 24 ]:

$\begin{array}{cc}\Omega \left(t\right)={\int }_0^{t}A\exp (-E_a/(\mathrm{R\prime T}\left(\mathrm{\tau \prime }\right)\left)\right)\mathrm{d\tau \prime }& \mathrm{(8)}\end{array}$

where A is the frequency factor, E_a is the activation energy, R′ is the universal gas constant (8.31 J/mol/K). T(τ) represents the absolute temperature history at past time τ′, which is the same physical quantity as the temperature field T in Eq. (5). Tissue is commonly considered necrotic when Ω≥1. For ex vivo porcine liver, we adopted literature-reported parameters in the HIFU temperature range (e.g., E_a=6.68×10⁵ J/mol, A=9.4×10¹⁰⁴ s⁻¹ [ 24 ]) and performed a sensitivity sweep (±20%) to assess robustness. We used the same simulated temperature–time fields to compute both CEM43 and Ω in parallel. We computed the thermal lesion using both the empirical CEM43 metric and the mechanistic Arrhenius integral. While CEM43 is convenient and widely used, the Arrhenius model is based on reaction kinetics and is better suited for predicting tissue damage under the high-temperature, short-duration conditions of HIFU. Therefore, we used the Arrhenius results to assess the validity of the CEM43 threshold (240 min) in our specific exposure scenario.

The acoustic characteristic parameters of biological tissues, such as attenuation coefficient, density, sound velocity, and nonlinear coefficient are related to temperature. When the temperature changes, the acoustic parameters will change accordingly, and the acoustic field distribution will also change accordingly. The change in the acoustic field will also cause the heat absorbed by the biological tissue to change, thereby affecting the temperature field in the biological tissue. Therefore, the acoustic field and temperature field generated by biological tissues under HIFU sonication will be coupled with each other. By dynamically calculating and analyzing the relationship between the harmonic ratio and the focal temperature and thermal lesion of porcine liver, an acoustic-thermal coupling model related to the harmonic ratio and thermal lesion is constructed, as shown in Figure 1. The acoustic field is solved using the finite difference frequency domain algorithm [ 18 , 25 ], and the temperature field is solved using the finite difference operator separation time domain algorithm [ 12 ].

Figure 1.Acoustic-thermal coupling model of harmonic ratio and thermal lesion. The model involves two main computational loops: the acoustic field is solved using the Finite-Difference Frequency-Domain (FDFD) method for the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation, and the temperature field is solved using the Finite-Difference Time-Domain (FDTD) method for the Pennes bioheat equation. The harmonic ratio `R<sup>*</sup>=(P<sub>3</sub>+P<sub>4</sub>)/P<sub>2</sub>` is calculated from the pressure amplitudes of the third (P₃) and fourth (P₄) harmonics, normalized by the second harmonic pressure P₂)

Numerical Implementation

1. Boundary Conditions

The acoustic field, governed by the KZK equation, utilized prescribed pressure P₀ at 1.4 MHz on the transducer surface (z=0), absorbing Perfectly Matched Layers (PML) at tissue boundaries to minimize reflections, and radial symmetry enforcement (∂p/∂r=0 at r=0). For the thermal field solved via Pennes equation, tissue boundaries maintained a fixed temperature of 30 °C with uniform initial conditions (T₀=30 ℃) throughout the domain.

2. Mesh Specifications

A structured grid covered the computational domain (axial span: 0–17 cm; radial span: 0–4 cm) with axial resolution Δz=0.1 mm and radial resolution Δr=0.05 mm, totaling 1,360,000 elements (1,700 axial × 800 radial nodes). The mesh distribution featured uniform spacing in the water region (0–11 cm) and graded refinement in the liver tissue (11–17 cm), with 2× higher nodal density near the focal zone (z=13 cm) to resolve nonlinear harmonics. Grid independence was confirmed through convergence validation, showing <1% variation in focal temperature when resolution was doubled (Δz′=0.05 mm, Δr′=0.025 mm). Convergence was ensured by progressively refining the mesh until the relative error in focal temperature fell below 1%. Specifically, the L2-norm error of the temperature field ‖T_h-T_h/2‖₂ /‖T_h/2‖₂ was computed between consecutive resolutions (where h denotes grid size), with convergence achieved at Δz=0.1 mm/ Δr=0.05 mm (error: 0.78%). The Courant–Friedrichs–Lewy (CFL) condition governed temporal stability, requiring acoustic time steps to satisfy Δt≤min(Δr,Δz)/c₀ (max CFL number: 0.89).

The HIFU irradiated biological tissue model is shown in Figure 2. The HIFU transducer is modeled after clinically approved concave spherical transducers (e.g., Chongqing Haifu JC200; Shenzhen Witasonic USgHIFU 1000), with a radius of a=5.5 cm, a geometric focal length of F=13 cm, and a center frequency of f=1.4 MHz. The transducer frequency (1.4 MHz) was selected to align with clinical HIFU standards for deep-tissue ablation [ 5 ], which optimizes the trade-off between penetration depth and focal precision in liver tissue. This frequency offers a balance between penetration depth and harmonic generation efficiency, as higher frequencies (e.g., 1.8 MHz) attenuate more rapidly, while lower frequencies (e.g., 1.0 MHz) generate weaker higher harmonics. At this frequency, attenuation in liver (α≈58 dB/m) allows sufficient energy deposition at the focal zone (depth: 11 cm water +3 cm tissue) while minimizing near-field heating. To assess the robustness of our findings with respect to this parameter, an additional simulation was conducted at 1.0 MHz and 1.8 MHz for sensitivity analysis, representing a frequency also commonly used in therapeutic ultrasound applications. The ultrasonic wave propagates Z_w=11 cm in water before reaching the tissue interface. Spatial peak pulse average intensity (I_SPPA) was used as the primary energy parameter, calculated from the acoustic pressure amplitude at the transducer surface (P₀) as I_SPPA= P₀²/(2ρ₀ c₀) [ 25 ]. To study nonlinear effects, I_SPPA ranged from 500 to 2000 W/cm² in increments of 250 W/cm². This range covers sub-boiling to boiling regimes relevant to clinical ablation [ 15 , 25 ]. The transducer and porcine liver tissue were placed in 30 ℃ water, and the porcine liver tissue dimensions were defined with L=6 cm along the acoustic axis (depth direction) and W=4 cm radially. Porcine liver was selected as the biological tissue model due to its well-established acoustic and thermal similarities to human liver tissue, including comparable attenuation coefficients, sound speeds, and nonlinearity parameters (values provided in Table 1) under HIFU exposure [ 24 , 26 , 27 ]. This computational model represents ex vivo homogeneous liver tissue, and thus the effects of blood perfusion and larger vasculature are not included. While this ex vivo simulation model provides a validated platform for studying fundamental acoustic-thermal interactions in homogeneous parenchymal tissue, we acknowledge inherent limitations in replicating in vivo physiological complexities (e.g., perfusion, microvasculature). Future studies will extend this framework to human-relevant in vivo models incorporating dynamic perfusion effects. This computational study focuses on fundamental thermal-acoustic interactions in homogeneous liver tissue, as explicitly stated in the description of the Pennes bioheat equation [ 20 ]: “This study focuses on homogeneous liver tissue without tumor-specific vasculature.” Consequently, tumor parameters are not incorporated here. The model geometry comprises two layers: (1) an 11 cm water layer (coupling medium) and (2) a 6 cm thick homogeneous liver tissue layer. Future tumor-specific models will incorporate multilayer architectures (e.g., skin/fat/muscle/tumor) with angiogenic parameters.

Figure 2. High-Intensity Focused Ultrasound (HIFU) sonicated a biological tissue model. Schematic diagram of the simulation setup showing the transducer and tissue geometry, where F is focal length, Z_w is water path length, L is tissue depth, W is tissue width, and a is transducer radius

The time steps of the acoustic field and temperature field are Δt_a=6.25×10^-8 s and Δt_t=60 ms, respectively. Spatial resolutions were selected to resolve the shortest harmonic wavelength (λ_min=c₀/(Nf₀)≈0.53 mm for N=4 at 1.4 MHz), with Δr=λ_min/10 ensuring ≥10 points per wavelength. The acoustic time step (Δt_a=6.25×10^-8 s) satisfied the CFL condition for numerical stability, while the thermal time step (Δt_t=0.06 s) was constrained by the tissue thermal relaxation time (τ~ρC_t Δz²/k_t≈ 1.3 s), ensuring Fourier number F₀<0.25) for accuracy [ 28 - 30 ]. Simulations were performed on a ThinkStation P620 (Windows 11 OS) with an Intel® Xeon® W-2295 CPU (18 cores, 3.0 GHz), 16 GB RAM, and NVIDIA RTX 5000 Graphics Processing Unit (GPU). While this computational study provides theoretical insights into harmonic ratio effects, future experimental validation using HIFU thermometry and tissue phantoms is planned to corroborate these findings. Table 1 gives the acoustic and thermal parameters of the relevant media.

3. Implementation Box

Implementation note for real-time monitoring is as follows: for translation to real-time systems, spectral windows of 5-10 ms with 50% overlap, Welch averaging (4-8 segments), and integration bands ±5% around 2f₀, 3f₀, and 4f₀ (where f₀ is the fundamental frequency) provide stable estimates of P₂, P₃, and P₄ at ≥100 Hz update rates. A coherence-weighted spatial average across a 10-15 mm aperture centered at the focus improves Signal-To-Noise Ratio (SNR) for R^*=(P₃+P₄)/P₂. The inferred temperature T=f(R^*, I, F₀) can be tabulated from our model across intensity (500-2000 W/cm²) and frequency (1.0-1.8 MHz) to enable sub-50 ms lookup during therapy.

Model Parameter Justification

The physical parameters for water and porcine liver used in this study were sourced from a comprehensive review of experimental literature to ensure the highest possible fidelity for an ex vivo model (Table 1). The temperature-dependent acoustic properties, particularly the attenuation coefficient critical for harmonic generation, were derived from polynomial fits to experimental data. These fits were based on extensive measurements of porcine liver tissue across a wide temperature range, as reported by Guntur et al. [ 31 , 32 ]. The use of this comprehensive, previously validated experimental dataset as the foundation for our model provides a high degree of confidence in its predictive accuracy, even in the absence of direct concurrent validation.

Dynamic Update of Temperature-Dependent Parameters

To capture the bidirectional acoustic-thermal coupling, the temperature-dependent acoustic attenuation coefficient of liver tissue, α(T), was dynamically updated at each computational time step based on the local instantaneous temperature. The relationship α(T) was implemented using a piecewise polynomial function derived from our previous experimental data fitting for ex vivo porcine liver [ 16 ]:

$\begin{array}{c}\alpha \left(T\right)=-\mathrm{2.217}\times {\mathrm{10}}^{\mathrm{-11}}\times T^7+\mathrm{9.868}\times {\mathrm{10}}^{\mathrm{-9}}\times T^6-\mathrm{1.842}\times {\mathrm{10}}^{\mathrm{-6}}\times T^5+\mathrm{0.0001844}\times T^4-\mathrm{0.01058}\times T^3+\mathrm{0.3481}\times T^2-\mathrm{6.189}\times T+\mathrm{50.7}\\ \mathrm{30\; °C}\le T\le \mathrm{90\; °C}& \mathrm{(9)}\end{array}$

for the temperature range 30 °C ≤T≤ 90 °C. This empirical function captures the characteristic decrease in attenuation from normal body temperature to approximately 50 °C, followed by a sharp increase above 50 °C due to protein denaturation and structural changes, as captured by the characteristic temperature-attenuation curve we previously reported in [ 16 ]). The values of the coefficients are: a₁=-2.217e^-11, a₂=9.868e^-9, a₃=-1.842e^-6, a₄=1.844e^-4, a₅=-0.01058, a₆=0.3481, a₇=-6.189, a₈=50.7. The nonlinearity parameter B/A was also considered temperature-dependent using a relation from [ 16 ], but its variation had a secondary effect compared to α(T) in the current model. This dynamic updating process is a key component of the coupling loop illustrated in Figure 1.

Results

By analyzing the relationship between harmonic ratio, focal temperature, and lesion size under varying HIFU sonication powers, the results highlight how nonlinear acoustic effects impact energy deposition and tissue damage.

Figure 3 shows the relationship between harmonic ratio and focal temperature at different HIFU sonication powers. The horizontal axis represents temperature and the vertical axis represents harmonic ratio. As shown in Figure 3, at different HIFU sonication powers, when the focal temperature is below 50 °C, the harmonic ratio gradually decreases with increasing temperature. When the focal temperature is above 50 °C, the harmonic ratio gradually increases with increasing temperature. While thermal denaturation in biological tissues generally initiates above approximately 42 °C, the specific turning point observed here at 50 °C pertains to the distinct temperature-dependent reversal of acoustic and thermal parameters in porcine liver tissue, particularly the attenuation coefficient, as evidenced by prior studies [ 33 , 34 ]. This 50 °C threshold marks the temperature where the rate of change of these parameters (e.g., the sign convention for how attenuation changes with temperature) reverses significantly in liver tissue [ 27 , 33 - 35 ], directly influencing harmonic generation efficiency and thus the observed turning behavior in the harmonic ratio. This phenomenon is specific to the tissue type and its characteristic parameter dependencies. While thermal denaturation in biological tissues generally initiates above approximately 42 °C, the specific turning point observed here at 50 °C pertains to the distinct temperature-dependent reversal of acoustic and thermal parameters in porcine liver tissue, particularly the attenuation coefficient, as evidenced by prior studies [ 34 - 37 ]. This 50 °C threshold marks the temperature where the rate of change of these parameters (e.g., the sign convention for how attenuation changes with temperature) reverses significantly in liver tissue [ 27 , 33 - 35 ], directly influencing harmonic generation efficiency and thus the observed turning behavior in the harmonic ratio. This reversal is primarily driven by the denaturation of collagen and other structural proteins, which alters the viscoelastic properties and hence the acoustic attenuation behavior of the tissue [ 34 , 36 ]. It is critical to note that the observed harmonic ratio-temperature relationship reflects a bidirectional acoustic-thermal coupling process. While harmonic generation enhances focal heating (acoustic energy conversion to thermal effects), rising temperature simultaneously alters tissue properties (e.g., attenuation coefficient), thereby modulating subsequent harmonic generation (thermal feedback to acoustic properties). The turning point at 50 °C emerges from this interplay, where thermal denaturation triggers a reversal in temperature-dependent acoustic parameters. This phenomenon arises from protein denaturation-induced structural reorganization in liver tissue. As temperatures approach 50 °C, the unfolding of collagen and other structural proteins [ 34 , 36 ] fundamentally alters the tissue microstructure. The disruption of hydrogen bonding and hydrophobic interactions modifies tissue viscoelasticity, thereby reversing the temperature dependence of key acoustic parameters (most notably the attenuation coefficient), as shown in references [ 34 , 36 ]. We acknowledge that inertial cavitation represents another important nonlinear phenomenon in HIFU therapy, known to influence energy deposition through bubble dynamics and shock wave generation [ 34 ]. However, this theoretical study specifically focuses on the role of harmonic ratio in thermal lesion formation via temperature-dependent acoustic parameter modulation. The exclusion of cavitation modeling is a deliberate scope limitation, as our coupled KZK-Pennes framework emphasizes thermally mediated harmonic generation mechanisms in homogeneous tissue. While cavitation effects were not computationally modeled here, their potential contribution to energy redistribution at high intensities (>1200 W/cm²) is noted in our interpretation of harmonic saturation effects [ 34 ]. Future work will incorporate cavitation models and experimental validation to quantify its relative contribution versus harmonic ratio effects. Regarding the harmonic intensity evolution, our simulations reveal that both second (P₂) and third (P₃) harmonic pressures exhibit progressive amplification with increasing HIFU intensity, adhering to the characteristic nonlinear scaling laws: P₂ demonstrates quadratic growth relative to fundamental intensity (P₂∝I²), while P₃ follows cubic scaling (P₃∝I³) in the pre-saturation regime. This intensification directly contributes to the harmonic ratio dynamics shown in Figure 3. However, beyond ~1200 W/cm² (corresponding to focal temperatures ≥50 °C), P₃ growth deviates from ideal cubic behavior due to energy redistribution to higher harmonics (n≥4) and temperature-mediated alterations in tissue nonlinearity and attenuation [ 34 , 37 ], consistent with established saturation phenomena in biological media. It is important to note that at these high intensities, the onset of inertial cavitation, which is not modeled in this study, would also contribute significantly to energy redistribution and absorption, further influencing the harmonic saturation effects observed. This complex interplay, particularly the transfer of energy from P₂ to P₃ and P₄, underscores the value of our (P₃+P₄)/P₂ ratio in capturing the dynamics of the thermally-dominant energy cascade. As reported by Guntur et al. in their experimental observations, at the turning temperature of 50 °C, the sign convention of the rate of change of thermal parameters with respect to temperature is reversed [ 31 , 32 ]. In addition, among all the acoustic and thermal parameters of porcine liver tissue, the attenuation coefficient is most significantly affected by temperature [ 35 ], and its amplitude turns around 50 °C [ 33 ].

Figure 3. The relationship between harmonic ratio value R^*=(P₃+P₄)/P₂ and focal temperature under different High-Intensity Focused Ultrasound (HIFU) sonication power

To elucidate this nonlinear energy transfer, the evolution of the second (P₂), third (P₃), and fourth (P₄) harmonic pressures was tracked and is presented in Figure 4d. The pressure amplitudes were normalized to their respective initial values (denoted as P_2,0, P_3,0, and P_4,0) to emphasize their relative growth dynamics. The vertical dashed line marks the time when the focal temperature reaches 50 °C, highlighting the turning point in the trends of the harmonic ratio and the attenuation coefficient. The initial decrease in the attenuation coefficient α from 30 °C to 50 °C (Figure 4c) reduces the absorption of higher harmonics, thereby moderating the rise of the harmonic ratio R^* (Figure 4b). Beyond 50 °C, the steep increase in α (Figure 4c) enhances the absorption of the higher harmonics (P₃ and P₄ in Figure 4d), leading to more efficient heating (Figure 4a) and a more efficient nonlinear energy cascade. This enhanced process is directly captured by the rapid increase in the harmonic ratio R^* (Figure 4b). This positive feedback loop is the fundamental mechanism behind the turning point.

Figure 4. Time-resolved evolution of (a) focal temperature, (b) harmonic ratio dynamics R^*=(P₃+P₄)/P₂, (c) dynamic attenuation coefficient α, and (d) normalized harmonic pressures P₂ (second harmonic), P₃ (third harmonic), and P₄ (fourth harmonic) during High-Intensity Focused Ultrasound (HIFU) sonication at I_SPPA=1500 W/cm². The vertical dashed line indicates the time when the focal temperature reaches 50 °C

Quantitatively, as the temperature increased from 50 ℃ to 70 ℃, the attenuation coefficient α increased by approximately 159.8% (from 44.3 dB/m to 115.1 dB/m, Eq. 8, Figure 4c). This directly led to a 7.5% increase in the harmonic ratio (Figure 4b and Figure 3) over the same temperature interval. This significant rise in attenuation is consistent with prior experimental observations of protein denaturation in liver tissue around 50 °C [ 31 , 32 ]. The relatively modest increase in the harmonic ratio likely reflects the competing effects of enhanced absorption of higher harmonics (which elevates the ratio) and potential saturation in nonlinear energy transfer at elevated temperatures. The growth rate of the third harmonic pressure (P₀) increased markedly after the turning point, as seen in Figure 4d, confirming the enhanced efficiency of the nonlinear energy transfer facilitated by the elevated attenuation.

Figure 5 shows the relationship between the harmonic ratio and the size of the lesion under different HIFU sonication powers. In Figure 5, the horizontal axis represents the harmonic ratio value, and the vertical axis represents the size of the lesion. As shown in Figure 5, as the harmonic ratio gradually increases, the thermal lesion also gradually increases. For lesion of the same size, the greater the HIFU sonication power, the smaller the value of the harmonic ratio, and for the same harmonic ratio value, the greater the HIFU sonication power, the greater the corresponding thermal lesion.

Figure 5. The relationship between harmonic ratio value R^*=(P₃+P₄)/P₂ and lesion size under different High-Intensity Focused Ultrasound (HIFU) sonication power

Frequency Sensitivity Analysis

To investigate the influence of the transducer frequency on the relationship between the harmonic ratio and focal temperature, we repeated the simulation using a 1.0 MHz and 1.8 MHz transducer, respectively, while keeping other parameters, such as the focal length and acoustic intensity, constant. Figure 6 compares the harmonic ratio versus focal temperature curves at 1.0 MHz, 1.4 MHz and 1.8 MHz for power of 105 W.

Figure 6. Sensitivity analysis of the harmonic ratio-temperature relationship R^*=(P₃+P₄)/P₂ at different transducer frequencies (1.0 MHz, 1.4 MHz, and 1.8 MHz) under power of 105 W

As shown in Figure 6, the characteristic turning-point behavior of the harmonic ratio at approximately 50 °C is robustly observed at both frequencies. This confirms that the phenomenon is a fundamental feature of the thermal-acoustic coupling process in liver tissue, rather than an artifact of a specific frequency. However, the frequency does modulate the characteristics of the curve. Specifically, the simulation at 1.8 MHz yields a higher overall harmonic ratio compared to 1.4 MHz and 1.0 MHz. This is expected, as higher frequencies lead to more efficient nonlinear energy transfer into higher harmonics and are also more readily absorbed, enhancing the effects captured by our defined ratio. This finding underscores the potential to tune the sensitivity of this biomarker by selecting the appropriate operating frequency. All runs used identical boundary conditions, mesh resolution, and power settings, ensuring that observed differences arise from frequency-dependent propagation and absorption. The persistence of the turning point across frequencies supports the robustness of the harmonic ratio as a biomarker, though its magnitude is frequency-dependent.

To address the impact of inertial cavitation on thermal lesion formation at high intensities (>1200 W/cm²), we conducted a comparative simulation incorporating a simplified cavitation model based on the Rayleigh–Plesset equation. This model accounts for bubble dynamics and associated heating effects, which are not included in the primary KZK-Pennes framework. The simulation was performed under the following conditions: I_SPPA=1500 W/cm², frequency f=1.4 MHz, sonication duration t=10 s, initial bubble radius R₀=2 μm, bubble density N_b=10¹⁰/m³, and tissue attenuation α=70.6 dB/m, consistent with the parameters used throughout the manuscript.

As illustrated in Figure 7, the inclusion of cavitation effects leads to a significantly larger and more irregular lesion boundary (red region, defined by the Arrhenius damage integral Ω≥1) compared to the lesion predicted by the purely thermal model (blue region). Quantitatively, cavitation increased the lesion cross-sectional area by approximately 68% (from 2.1 mm² to 3.5 mm²) and reduced the time to reach the critical lesion threshold (Ω=1) from 4.2 s to 2.8 s. This accelerated and enlarged lesion formation is attributable to the additional heat deposition from bubble collapse, microstreaming, and shockwave emission mechanisms not captured by the thermal model alone. These results underscore the substantial influence of cavitation on energy deposition and lesion morphology under high-intensity conditions. While the primary focus of this study is on the harmonic ratio as a biomarker for thermally mediated nonlinear effects, this supplementary analysis confirms that cavitation becomes a dominant factor at intensities above 1200 W/cm².

Figure 7. Comparative simulation of thermal lesion formation with and without cavitation effects at an intensity of 1500 W/cm² and 10 s sonication duration. The lesion boundary is defined by the Arrhenius damage integral Ω≥1

In response to the applicability of the Sapareto–Dewey thermal dose formulation under rapid heating conditions, we performed a comparative analysis of lesion boundaries predicted by both the empirical Sapareto–Dewey model (CEM43≥240 min) and the mechanistic Arrhenius integral model (Ω≥1). The simulation was conducted at a representative intensity of 1000 W/cm² with a 6 s sonication duration under ex vivo conditions (no perfusion, no cavitation), ensuring monotonic heating conducive to both models.

To further evaluate the consistency between empirical and mechanistic thermal dose models, we performed a comparative analysis of lesion boundaries predicted by both the empirical Sapareto–Dewey model (based on thermal dose CEM43≥240 min) and the mechanistic Arrhenius integral model (based on Ω≥1). The simulation was conducted at a representative intensity of 1000 W/cm² with a 6 s sonication duration under ex vivo conditions (no perfusion, no cavitation), ensuring monotonic heating conducive to both models. This intensity level was chosen to induce substantial thermal lesioning while remaining below the threshold for significant boiling or cavitation, allowing a clear interrogation of model agreement under controlled, sub-boiling conditions.

Figure 8 illustrates this comparison. As shown in Figure 8, the lesion boundaries predicted by both models exhibit strong agreement. Subplot (a) shows the lesion defined by CEM43≥240 min, subplot (b) shows the lesion defined by Ω≥1, and subplot (c) illustrates their high degree of overlap. Quantitative analysis revealed that the Arrhenius-based lesion volume differed from the CEM43-based volume by only 8.5%, with a Dice similarity coefficient of 0.94. This close agreement validates the use of the CEM43 metric under the controlled, sub-boiling conditions simulated here, where heating is rapid yet monotonic.

Figure 8.Comparative analysis of thermal lesion boundaries predicted by Sapareto-Dewey and Arrhenius models at 1000 W/cm². (a) Lesion defined by thermal dose CEM43≥240 min, (b) lesion defined by Arrhenius damage integral Ω≥1, (c) overlay of the two lesion boundaries (CEM43: Cumulative Equivalent Minutes at 43 °C)

It is noteworthy that while both models converge well under these conditions, the kinetically-based Arrhenius model provides a more fundamental description of tissue lesion and is considered more appropriate for predicting lesion formation under extreme exposure protocols characterized by very rapid heating rates (>50 °C/s) or peak temperatures exceeding 80 °C, where the empirical CEM43≥240 min formulation may be extrapolated beyond its validated range. Therefore, while CEM43≥240 min remains a practical and widely adopted metric in HIFU literature, we recommend the use of the Arrhenius integral for scenarios involving complex thermal histories or very short-duration, high-temperature exposures.

Discussion

The findings of this study provide significant insights into the influence of harmonic ratio on HIFU-induced thermal lesion formation in biological tissue, offering a deeper understanding of the nonlinear acoustic effects that govern HIFU therapy. The results demonstrate that harmonic ratio plays a critical role in the spatial distribution of acoustic energy deposition and the resulting thermal lesion size, which is consistent with previous research on the nonlinear behavior of ultrasound in biological tissues [ 34 , 35 ]. Our results demonstrate that the influence of the R^*=(P₃+P₄)/P₂ ratio extends beyond previous findings related to simpler ratios like P₂/P₁ [ 16 ], by directly linking the dynamics of higher-order harmonics (P₃ and P₄) to the most thermally significant part of the energy deposition process, as evidenced by the identified turning point at 50 °C. As shown in Figure 3, the relationship between harmonic ratio and focal temperature exhibits a distinct turning point at 50 °C, which corresponds to the thermal denaturation threshold of liver tissue. The turning point in harmonic ratio behavior at 50 ℃ is attributed to the onset of thermal denaturation in porcine liver tissue. Our model, through dynamic updating of the temperature-dependent attenuation coefficient α(T) (Eq. 8, Figure 4c), explicitly captures this process. This simulated attenuation trend, which reverses (decreasing then increasing sharply above 50 °C) in agreement with experiments [ 32 , 34 , 36 ], constitutes the fundamental reason why the harmonic ratio exhibits a parallel reversal. Our time-resolved simulations (Figure 4) numerically demonstrate this bidirectional acoustic-thermal coupling: the changing attenuation coefficient α(T) directly alters the efficiency of harmonic generation and absorption (Figure 4d), which in turn modulates the localized heating rate (Figure 4a). Therefore, the 50 ℃ turning point is not merely an empirical observation but a direct consequence of a fundamental, denaturation-induced change in the tissue’s acoustic properties. This phenomenon is primarily driven by protein denaturation (notably collagen unfolding) around this temperature [ 34 , 36 ], which fundamentally reorganizes the tissue microstructure. This structural alteration disrupts intermolecular forces (e.g., hydrogen bonding, hydrophobic interactions), leading to a reversal in the temperature dependence of viscoelastic properties and, consequently, key acoustic parameters such as the attenuation coefficient [ 33 , 34 , 36 , 37 ]. This bidirectional energy transfer (acoustic energy conversion to heat followed by heat-induced alteration of acoustic properties) explains why harmonic ratio both influences and is influenced by temperature, reconciling apparent contradictions in causal directionality. The efficacy of the harmonic ratio as a sensitive indicator of this shift is underscored by its direct reflection of the nonlinear energy cascade efficiency into higher harmonics. This cascade is governed by the attenuation coefficient and the nonlinearity parameter B/A, both of which exhibit strong temperature dependence [ 33 , 34 ]. The initial decrease in attenuation up to ~50 °C reduces the absorption of higher harmonics, thereby tempering the rise of the harmonic ratio. Beyond this critical temperature, protein denaturation induces a substantial increase in attenuation, which preferentially enhances the absorption of higher-order harmonics (e.g., P₃, P₄), consequently amplifying the harmonic ratio. This creates a positive feedback loop that intensifies focal heating. Thus, the 50 °C turning point, captured by our dynamically coupled model, is mechanistically explained by thermally induced structural reorganization that fundamentally alters the tissue’s acoustic properties, validating the ratio’s role as a biomarker. The observed behavior of harmonic ratio at different HIFU sonication powers suggests its potential as a real-time temperature monitoring tool during HIFU treatment. This capability could enhance the precision of thermal control, ensuring that the target tissue reaches the desired temperature while minimizing lesion to surrounding healthy tissue.

Practical framework for real-time harmonic-ratio–guided feedback in clinical HIFU

To translate the harmonic ratio biomarker into clinical practice, we propose a closed-loop control architecture wherein R^*=(P₃+P₄)/P₂ serves as a real-time surrogate for the focal thermal state. The core of this system is a Passive Acoustic Monitoring (PAM) array, integrated within the HIFU transducer housing, featuring 64 to 256 elements with a 1–5 MHz bandwidth and a pitch of 0.3–0.5 mm. This array acquires backscattered and radiated acoustic emissions during sonication. Each channel is sampled at 20–40 MS/s, facilitated by low-noise preamplifiers and analog anti-aliasing filters (cutoff ≈8–10 MHz) to ensure a sufficient signal-to-noise ratio for resolving the second through fourth harmonics.

Real-time processing is implemented on an onboard GPU or Field-Programmable Gate Array (FPGA), which performs Short-Time Fourier Transform (STFT) using 5–10 ms windows with 50% overlap, yielding an update rate of 100–200 Hz. Spectral variance is reduced using a Welch or multi-taper method, and the pressure amplitudes P₂, P₃, and P₄ are obtained by integrating power within ±5% of each harmonic’s center frequency, with outliers rejected via a Hampel filter. The computed harmonic ratio R^* is spatially aggregated across a focal aperture using coherence weighting. Focal temperature is inferred through a precomputed lookup table (i.e., a data structure for fast temperature estimation) T=f(R^*, I, F₀) derived from our coupled KZK–Pennes model or, alternatively, via an online extended Kalman filter that fuses R^* with a simplified thermal state model. A Model Predictive Controller (MPC) subsequently modulates the acoustic intensity and duty cycle to track a predefined target temperature trajectory (e.g., a 60–65 °C plateau). This control loop operates under stringent safety constraints, including upper thresholds for R^* (indicative of rapid heating), cumulative thermal dose (CEM43 or Arrhenius integral Ω), and cavitation monitors (e.g., a rise of more than 6 dB in broadband high-frequency energy), triggering immediate power reduction or cessation if violated.

The feasibility of this framework is supported by practical computational estimates; a system with 128 channels sampling at 20 MS/s requires approximately 2–3 GFLOPs (Giga Floating-Point Operations Per Second) per frame for Welch-based spectral averaging, which a high-performance GPU (e.g., NVIDIA RTX series) can process at over 200 frames/s, resulting in an end-to-end latency of <50 ms. For hepatic targets at 10–13 cm depth, published array sensitivities suggest PAM SNRs exceeding 12 dB for P₂ and 8–10 dB for P₃/P₄ at intensities of 1000–1500 W/cm², permitting stable estimation of R^* with <10% variance after temporal smoothing. Integration with clinical HIFU systems would proceed through an auxiliary PAM front-end, with validation spanning tissue-mimicking phantoms, ex vivo liver studies, and ultimately in vivo large-animal models to benchmark tracking accuracy, control performance, and safety. Finally, regulatory and cybersecurity imperatives necessitate that the feedback module operates as an isolated decision-support layer with auditable logs, deterministic fail-safes, and cryptographically signed software integrity, ensuring robustness within the clinical environment.

The proposed framework will be rigorously validated through benchtop experiments, ex vivo tissue studies, and eventually in vivo trials. Key validation metrics will include temperature tracking Root Mean Square Error (RMSE) against thermometry (<1.5 °C target), control overshoot (<2 °C), settling time (<3 s), and safety sentinel false-negative rate (<1%). These metrics ensure that the harmonic-ratio-guided feedback system meets clinical safety and efficacy standards.

The correlation between harmonic ratio and lesion size, as illustrated in Figure 5, further underscores the importance of harmonic ratio in HIFU therapy. Higher harmonic ratio values were associated with larger thermal lesion, indicating that harmonic ratio components enhance energy deposition in the focal region, leading to more extensive tissue lesion. This strong positive correlation validates our hypothesis that the R^*=(P₃+P₄)/P₂ ratio serves as an effective proxy for the efficiency of localized, thermally-significant energy deposition. This finding aligns with the hypothesis that nonlinear acoustic effects, including harmonic ratio generation, play a crucial role in the efficiency of HIFU-induced thermal ablation. Notably, an inverse correlation between HIFU power and harmonic ratio values emerges when comparing lesions of identical size. This suggests that at elevated power levels, the efficiency of the specific harmonic cascade captured by our ratio may decrease. This could be due to the saturation of the harmonic generation process itself, or more likely, the activation of competing and highly efficient energy dissipation mechanisms. Chief among these is inertial cavitation, which becomes a dominant factor at high intensities. Cavitation bubbles can absorb and scatter acoustic energy, and their violent collapse generates localized shock waves and microstreaming, creating an additional, potent heat source not accounted for in our current thermal-harmonic model [ 34 , 37 ]. At intensities exceeding 1200 W/cm², inertial cavitation likely becomes a significant heat source, potentially altering the harmonic ratio behavior. This suggests that higher sonication power levels may suppress harmonic ratio generation, possibly due to increased nonlinear effects or cavitation, which can alter the acoustic field and energy distribution [ 34 , 37 ]. While reference [ 37 ] valuably characterized the intensity threshold (approximately 8 W/cm² at 4 MHz) for significant nonlinear waveform distortion in multilayer superficial tissues (skin/fat/muscle), our study diverges fundamentally in its core objective, biological context, and key findings. We focus on establishing the harmonic ratio (P₃+P₄)/P₂ as a quantitative biomarker and elucidating its dynamic, bidirectional coupling with temperature specifically identifying a distinct turning point at 50 °C linked to thermal denaturation in homogeneous liver tissue (Figure 3). Our study simulates this temperature-modulated harmonic behavior using the KZK equation, which is optimized for deep-tissue ablation at 1.4 MHz in the liver. This approach differs from that of reference [ 37 ], which employed the Westervelt equation to analyze pressure and temperature rises due to the onset of nonlinearity in adipose tissue, focusing on superficial applications. Consequently, the primary therapeutic parameter we advance is the harmonic ratio itself. We correlate this ratio with lesion size (Figure 5) to establish its potential for real-time monitoring in hepatic HIFU therapy, thereby complementing reference [ 37 ]’s focus on intensity thresholds for superficial treatments. This observation highlights the need for careful optimization of HIFU parameters, such as power and exposure time, to achieve the desired therapeutic outcome while minimizing unintended tissue lesion.

Furthermore, our sensitivity analysis demonstrated that while the fundamental turning-point behavior of the harmonic ratio is robust across different frequencies (Figure 6), the choice of frequency significantly modulates the magnitude of the harmonic ratio. The observed increase in the harmonic ratio at 1.4 MHz compared to 1.0 MHz and 1.8 MHz aligns with nonlinear acoustic theory, which predicts more pronounced harmonic generation at higher frequencies for a given pressure amplitude. This frequency-dependent behavior is clinically relevant, as it suggests that the harmonic ratio could be optimized as a monitoring biomarker by selecting a frequency that provides the best signal-to-noise ratio for detecting thermal changes. It also highlights the complexity of HIFU treatment planning, where frequency must be balanced to achieve sufficient penetration depth, precise focusing, and efficient nonlinear heating.

When considering feasibility and robustness, the harmonic ratio R^*=(P₃+P₄)/P₂ is inherently a spectral-derived metric. Its value is determined by the relative energy distribution between harmonics, making it robust to common confounding factors that affect absolute signal amplitudes, such as transducer output drift or tissue transmission losses. Residual risks include confounding from cavitation-induced broadband energy; this is mitigated by concurrent cavitation sentinels and by excluding broadband bins from the harmonic integrals. In low-SNR scenarios (e.g., deep targets or high adiposity), temporal smoothing (100-200 ms), aperture compounding, and adaptive gain control stabilize R^*=(P₃+P₄)/P₂ without materially increasing latency. If any quality metric (e.g., SNR <6 dB, coherence <0.5) deteriorates, the controller reverts to a conservative open-loop plan with reduced power until signal quality recovers.

Regarding interpretability and clinical readiness, it might be tempting to infer that any rise in harmonic content directly guarantees sufficient thermal dose, but the mapping between spectral changes and tissue injury depends on intensity, frequency, and tissue pathway. Our framework addresses this by calibrating R^*-to-temperature mappings under the specific therapy protocol, enforcing independent safety sentinels, and reverting to conservative operation when signal quality degrades. This layered design aligns with best practices for safety-critical ultrasound systems and avoids over-reliance on a single biomarker.

The potential of harmonic ratio as a biomarker for assessing tissue lesion during HIFU therapy is another important implication of this study. By monitoring harmonic ratio levels, clinicians could evaluate the extent of tissue lesion in real-time and adjust treatment parameters accordingly. This capability is particularly relevant in clinical settings where precision is critical, such as in the treatment of liver or prostate tumors. This aligns with findings from studies using passive acoustic mapping, which have demonstrated that changes in harmonic content can be correlated with the spatial extent of tissue necrosis. The ability to assess tissue lesion non-invasively using harmonic ratio could significantly improve the safety and efficacy of HIFU therapy, reducing the risk of collateral lesion to surrounding healthy tissue.

The results of this study also have broader implications for the optimization of HIFU therapy. The observed relationship between harmonic ratio and thermal lesion formation suggests that controlling harmonic ratio generation could be a key strategy for enhancing treatment outcomes. For example, by modulating HIFU parameters to promote or suppress harmonic ratio generation, clinicians could tailor the treatment to achieve the desired balance between therapeutic efficacy and safety. This approach could be particularly beneficial in cases where precise control of the thermal lesion size is critical, such as in the treatment of small or irregularly shaped tumors.

This study, while illuminating the fundamental relationship between harmonic ratio and thermal lesion formation, has several limitations that pave the way for future research. First, the model’s physiological relevance is limited by the exclusion of blood perfusion in the Pennes bioheat equation. This simplification omits the significant cooling effect of blood flow, which would alter the absolute temperature values and lesion growth rates in vivo, particularly in highly vascularized tissues or tumors. Second, and critically, this study deliberately excludes the effects of inertial cavitation to isolate the thermal-acoustic feedback loop mediated by the harmonic ratio. At the higher intensities investigated (>1200 W/cm²), cavitation is known to be a major contributor to both heating and mechanical tissue damage. A more comprehensive model would require coupling the current framework with a bubble dynamics model, such as the Rayleigh-Plesset equation, to account for the additional heat source from cavitation activity. To provide a quantitative estimate of this effect’s potential magnitude, we conducted a preliminary simulation comparing lesion formation with and without a simplified cavitation heat source. As shown in Figure 7, the inclusion of cavitation can significantly increase the predicted lesion volume, underscoring the importance of incorporating this mechanism in future, more comprehensive models. Such an extension would allow for a direct comparison of lesion formation with and without cavitation effects, which is a primary goal of our future work. This will enable a more accurate prediction of total lesion volume and shape under clinically relevant high-power sonication protocols. Third, the biological environment is simplified by considering homogeneous liver tissue, neglecting the complex layered architecture (e.g., skin, fat, muscle) and the specific angiogenic properties of tumors. Fourth, the findings are based on a purely computational model and await experimental validation. Finally, the immediate next step is to conduct dedicated ex vivo experimental validation using a calibrated HIFU system to crucially test the predictive power of the harmonic ratio biomarker.

These limitations define clear directions for future work. Firstly, to address perfusion, the model will be extended by incorporating the full Pennes perfusion term (ω_b, ρ_b, c_b,(T_α-T)). Sensitivity analyses will be conducted across a range of perfusion rates (e.g., 0-20 kg/m³/s) representative of normal liver parenchyma and tumor tissue to quantitatively assess its impact on the harmonic ratio-thermal lesion relationship. Secondly, the model will evolve to include multi-layered anatomical structures with appropriate acoustic and thermal properties for each layer to enhance clinical translatability for specific applications (e.g., abdominal or breast tumor ablation). Fourth, this study focused on establishing the fundamental behavior of the harmonic ratio at representative frequencies. A comprehensive parametric sweep, including a wider range of frequencies and intensities, would be necessary to build a complete optimization framework for HIFU therapy based on this biomarker. Such a study would aim to identify optimal parameter ranges for specific clinical scenarios, which represents a significant and valuable direction for future research. Applicability of thermal lesion models: In our homogeneous, ex vivo liver simulations with monotonic temperature rise, CEM43 and Arrhenius criteria produced closely matching lesion predictions. However, in scenarios with highly transient heating such as boiling onset, cavitation-enhanced energy deposition, rapid cooling from perfusion, or multi-pulse duty cycles with recovery, Arrhenius kinetics should be the primary criterion, with CEM43 reported as a secondary, interpretability-focused metric. Future in vivo studies will therefore adopt Arrhenius as the default lesion model and use CEM43 for clinical comparability, accompanied by uncertainty quantification and parameter sensitivity analysis (A, E_a). Clarification on scope and monitoring roles: While the harmonic ratio provides a sensitive thermal-state surrogate, it is not intended to replace thermometry or cavitation monitoring in complex in vivo scenarios. Rather, it complements these modalities within a layered safety architecture. In practice, R^*-based control should operate alongside thermometry (e.g., Proton Resonance Frequency Shift (PRFS) Magnetic Resonance Imaging (MRI) or ultrasonic thermometry when available) and cavitation monitoring, with predefined arbitration rules that prioritize patient safety. Finally, and most importantly, the immediate next step is to conduct dedicated ex vivo experimental validation using a calibrated HIFU system. This will involve irradiating fresh porcine liver samples while simultaneously monitoring focal temperature (e.g., with thermocouples or MRI thermometry) and acoustic emissions (with a hydrophone) to measure the harmonic ratio empirically. Subsequent histopathological analysis (e.g., Triphenyl Tetrazolium Chloride (TTC) of the induced lesions will allow for a direct, quantitative comparison between the simulated and measured outcomes, crucially testing the predictive power of the harmonic ratio biomarker.

Conclusion

This study explored the influence of a novel higher-order harmonic ratio R^*=(P₃+P₄)/P₂ on HIFU-induced thermal lesion formation in biological tissue using a coupled acoustic-thermal model. The use of this ratio, combined with a dynamically coupled KZK-Pennes model, provides a more sensitive and mechanistic understanding of the acoustic-thermal feedback process than previously achieved with second-harmonic metrics alone. The results revealed that harmonic ratio significantly impacts the spatial distribution of acoustic energy deposition and the resulting thermal lesion size. Higher harmonic ratio values correlated with larger lesion, while lower values led to smaller lesion, highlighting the importance of controlling harmonic ratio generation to optimize HIFU therapy. A critical turning point was observed at 50 °C, where harmonic ratio behavior changed due to thermal denaturation of liver tissue, suggesting its potential as a sensitive indicator for real-time temperature monitoring during treatment. Furthermore, the correlation between harmonic ratio and lesion size indicates its utility as a biomarker for assessing tissue lesion. These findings underscore the potential of harmonic ratio in enhancing the precision and safety of HIFU therapy. Future research should focus on elucidating the mechanisms of harmonic ratio generation and its relationship with HIFU parameters, such as frequency and exposure time, to further refine treatment strategies and improve therapeutic outcomes.

Acknowledgment

The authors sincerely thank the anonymous reviewers for their helpful comments and suggestions.

Authors’ Contribution

D. Hu conceived and designed the research idea and framework, and wrote the manuscript; H. Jiwen performed the simulations; Ch. Wei revised the paper. All authors have read and agreed to the published version of the manuscript.

Ethical Approval

The School of Information Science and Engineering at Changsha Normal University has approved the research and use of necessary resources.

Funding

This study was funded by the Hunan Provincial Department of Education Key Project (25A0732), and Changsha Natural Science Foundation Project (grant No.kq2202313).

Conflict of Interest

None

References

1. Ashar H, Ranjan A. Immunomodulation and targeted drug delivery with high intensity focused ultrasound (HIFU): Principles and mechanisms. Pharmacol Ther. 2023; 244:108393. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
2. Sofuni A, Takeuchi H, Sugimoto K, Itoi T, Miyazawa H. High-intensity focused ultrasound treatment for hepatocellular carcinoma. J Med Ultrason (2001). 2024;1-5. DOI | PubMed
3. Calik J, Zawada T, Bove T, Dzięgiel P, Pogorzelska-Antkowiak A, Mackiewicz J, Woźniak B, Sauer N. Healing process after high-intensity focused ultrasound treatment of benign skin lesions: dermoscopic analysis and treatment guidelines. J Clin Med. 2024; 13(4):931. DOI
4. Zulkifli D, Manan HA, Yahya N, Hamid HA. The Applications of High-Intensity Focused Ultrasound (HIFU) Ablative Therapy in the Treatment of Primary Breast Cancer: A Systematic Review. Diagnostics (Basel). 2023; 13(15):2595. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
5. Zhou K, Strunk H, Dimitrov D, Vidal-Jove J, Gonzalez-Carmona MA, et al. US-guided high-intensity focused ultrasound in pancreatic cancer treatment: a consensus initiative between Chinese and European HIFU centers. Int J Hyperthermia. 2024; 41(1):2295812. DOI | PubMed
6. Roohi R, Baroumand S, Hosseinie R, Ahmadi G. Numerical simulation of HIFU with dual transducers: The implementation of dual-phase lag bioheat and non-linear Westervelt equations. Int J Heat Mass Transf. 2021; 120:105002. DOI
7. Rahpeima R, Lin CA. A comprehensive numerical procedure for high-intensity focused ultrasound ablation of breast tumour on an anatomically realistic breast phantom. PLoS One. 2024; 19(10):e0310899. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
8. Yadav SK, Paul S, Singh MS. Effect of HIFU-Induced Thermal Ablation in Numerical Breast Phantom. Photonics. 2023; 10(4):425. DOI
9. Chen J, LeBlang S, Hananel A, Aginsky R, Perez J, Gofeld M, et al. An incoherent HIFU transducer for treatment of the medial branch nerve: Numerical study and in vivo validation. Int J Hyperthermia. 2020; 37(1):1219-28. DOI | PubMed
10. Huang W, Jiao Y, Li J, He Y, Shao W, Cui Y. Evaluation of Dual-Frequency Switching HIFU for Optimizing Superficial Ablation. Ultrasound Med Biol. 2024; 50(6):908-19. DOI | PubMed
11. Zou X, Qian S, Tan Q, Dong H. Formation of thermal lesions in tissue and its optimal control during HIFU scanning therapy. Symmetry. 2020; 12(9):1386. DOI
12. Kim SJ, Hwang JY, Kim YJ, Pae KN. Numerical simulation method for prediction of HIFU induced lesions in human tissue: FDTD-LBM. Phys Wave Phen. 2023; 31(1):30-5. DOI
13. Fura L, Zolek N, Kujawska T. Numerical simulations and experimental verification of the extent of HIFU-induced tissue necrosis. In Medical Imaging 2022: Image-Guided Procedures, Robotic Interventions, and Modeling; San Diego, California, United States: SPIE; 2022. p. 634-8.
14. Fura L, Tymkiewicz R, Kujawska T. Numerical studies on shortening the duration of HIFU ablation therapy and their experimental validation. Ultrasonics. 2024; 142:107371. DOI | PubMed
15. Gharloghi S, Gholami M, Haghparast A, Dehlaghi V. Numerical study for optimizing parameters of high-intensity focused ultrasound-induced thermal field during liver tumor ablation: HIFU Simulator. Iran J Med Phys. 2017; 14(1):15-22. DOI
16. Dong H, Liu G, Tong X. Influence of temperature-dependent acoustic and thermal parameters and nonlinear harmonics on the prediction of thermal lesion under HIFU ablation. Math Biosci Eng. 2021; 18(2):1340-51. DOI | PubMed
17. Zhou Y, Yu Z, Ma Q, Guo G, Tu J, Zhang D. Noninvasive Treatment-Efficacy Evaluation for HIFU Therapy Based on Magneto-Acousto-Electrical Tomography. IEEE Trans Biomed Eng. 2019; 66(3):666-74. DOI | PubMed
18. Li C, Chen S, Wang Q, Li H, Xiao S, Li F. Effects of thermal relaxation on temperature elevation in ex vivo tissues during high intensity focused ultrasound. IEEE Access. 2020; 8:212013-21. DOI
19. Richards N, Christensen D, Hillyard J, Kline M, Johnson S, Odéen H, Payne A. Evaluation of acoustic-thermal simulations of in vivo magnetic resonance guided focused ultrasound ablative therapy. Int J Hyperthermia. 2024; 41(1):2301489. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
20. Suarez-Castellanos IM, De Sallmard G, Vanstaevel G, Ganeau A, Bawiec C, Chapelon JY, et al. Dynamic Ultrasound Focusing and Centimeter-Scale Ex Vivo Tissue Ablations With a CMUT Probe Developed for Endocavitary HIFU Therapies. IEEE Trans Ultrason Ferroelectr Freq Control. 2023; 70(11):1470-81. DOI | PubMed
21. Wessapan T, Rattanadecho P. Flow and heat transfer through a porous tumor during high-intensity focused ultrasound. Int J Heat Mass Transf. 2023; 216:124501. DOI
22. Lari S, Han SW, Kim JU, Kwon HJ. Design of HIFU treatment plans using thermodynamic equilibrium algorithm. Algorithms. 2022; 15(11):399. DOI
23. Padilla F, Ter Haar G. Recommendations for Reporting Therapeutic Ultrasound Treatment Parameters. Ultrasound Med Biol. 2022; 48(7):1299-308. DOI | PubMed
24. Tan Q, Zou X, Ding Y, Zhao X, Qian S. The influence of dynamic tissue properties on HIFU hyperthermia: A numerical simulation study. Appl Sci. 2018; 8(10):1933. DOI
25. Gu J, Jing Y. Modeling of wave propagation for medical ultrasound: a review. IEEE Trans Ultrason Ferroelectr Freq Control. 2015; 62(11):1979-93. DOI | PubMed
26. Bawiec CR, Rosnitskiy PB, Peek AT, Maxwell AD, Kreider W, Haar GRT, et al. Inertial Cavitation Behaviors Induced by Nonlinear Focused Ultrasound Pulses. IEEE Trans Ultrason Ferroelectr Freq Control. 2021; 68(9):2884-95. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
27. Civale J, Rivens I, Shaw A, Ter Haar G. Focused ultrasound transducer spatial peak intensity estimation: a comparison of methods. Phys Med Biol. 2018; 63(5):055015. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
28. LeVeque RJ. Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems. Society for Industrial and Applied Mathematics; 2007.
29. Polycarpou AC. Introduction to the finite element method in electromagnetics. Springer; 2006.
30. Clarke RL, Bush NL, Ter Haar GR. The changes in acoustic attenuation due to in vitro heating. Ultrasound Med Biol. 2003; 29(1):127-35. DOI | PubMed
31. Guntur SR, Choi MJ. Influence of temperature-dependent thermal parameters on temperature elevation of tissue exposed to high-intensity focused ultrasound: numerical simulation. Ultrasound Med Biol. 2015; 41(3):806-13. DOI | PubMed
32. Guntur SR, Lee KI, Paeng DG, Coleman AJ, Choi MJ. Temperature-dependent thermal properties of ex vivo liver undergoing thermal ablation. Ultrasound Med Biol. 2013; 39(10):1771-84. DOI | PubMed
33. Dong H, Xiao Z, Qian S. Simulation study on the influence of temperature-dependent acoustic-thermal parameters on tissue lesion under HIFU irradiation. UPB Sci Bull Series A. 2020; 82(2):207-20.
34. Xu P, Wu H, Shen G. Characterization of weakly nonlinear effects in relationship to transducer parameters in focused ultrasound therapy. Med Phys. 2024; 51(10):7619-731. DOI | PubMed
35. Han M, Song W, Zhang F, Li Z. Modeling for Quantitative Analysis of Nakagami Imaging in Accurate Detection and Monitoring of Therapeutic Lesions by High-Intensity Focused Ultrasound. Ultrasound Med Biol. 2023; 49(7):1575-85. DOI | PubMed
36. Sanderson SG, Easthope B, Farias C, Doddridge I, Cook JA, Dahl DB, Dillon CR. Characterizing temperature-dependent acoustic and thermal tissue properties for high-intensity focused ultrasound computational modeling. Int J Thermophys. 2024; 45(10):143. DOI
37. Mortazavi SMJ, Mokhtari-Dizaji M. Threshold of Linear and Non-Linear Behavior of High Intensity Focused Ultrasound (HIFU) in Skin, Fat, and Muscle Tissue Using Computer Simulation. Iran J Med Phys. 2022; 19(3):181-8. DOI