Numerical evaluation of high-intensity focused ultrasound- induced thermal lesions in atherosclerotic plaques

1.   Introduction

Atherosclerosis is a chronic, inflammatory disease and starts as an endothelial dysfunction, which initiates recruitment of inflammatory cells and lipids within the arterial wall ^[1]. The main species and substances play an important role in early atherosclerosis development, i.e., LDL, oxidized LDL, monocytes, macrophages, foam cells, smooth muscle cells, cytokines and collagen ^[2]. During the progression of atherosclerotic plaque, macrophages play a critical role in pro-inflammatory regulation and lipoprotein homeostasis ^[3]. This suggests that macrophage denaturation and ablation by hyperthermia may actively suppress inflammation through inhibiting the release of pro-inflammatory cytokines from macrophages and effectively prevent the recruitment of macrophages in plaque ^[3,4]. On the other hand, collagen fibrin is the end product of blood clotting/atherosclerotic plaque that constitutes a proteinaceous 3D network providing a filamentous mechanical scaffold of clots and thrombi ^[5]. Denaturation of proteins in the membrane will inhibit of the activity of enzymes ^[6].

Generally, Macromolecular and cell will start to denature between 44 ℃ ^[7] and 47 ℃ ^[8,9]. If the temperature rise is above a threshold of 55 ℃ and the exposure time is at least 1 second, irreversible cell death will be induced through coagulation ablation ^[10,11]. Meanwhile, the equilibrium of chemical concentration in atherosclerotic plaque is destroyed and a high reduction in enzyme activity is achieved ^[12,13]. Therefore, for temperatures above 60 ℃, biochemical reactions regulate the oxidation of LDL, involving enzymes (such as Lp-LpA2, cysteine-, serine- and metallopeptidases) and free radicals in the endothelium greatly reduced and cause macrophage death and reduced enzymes activity in atherosclerosis ablation ^[14]. Coagulation necrosis of macrophages may stop the ingestion of large amounts of ox-LDL by the fatty macrophages for transforming into foam cells, and thus slow the growth of atherosclerotic plaque ^[14,15].

High-intensity focused ultrasound (HIFU) has been shown to produce coagulative necrosis and is currently used for clinical medical without significant biological damage to intervening and surrounding tissue ^[16,17]. Compared to the traditional approaches for atherosclerosis treatment such as Balloon Angioplasty ^[18], atherectomy ^[19], and surgical bypass ^[20,21] etc., the major advantages of HIFU therapy are its noninvasive and nonionizing characteristics, as well as the absence of long-term cumulative effects, and fast, confined lesion formation ^[22,23], and can be controlled selectively to obtain a desired temperature distribution and elevation in the treatment area with proper choice of operating frequency and transducer design ^[24,25]. Siegel et al. performed one of the first feasibility studies of HIFU for ablation of atherosclerotic plaques in human atherosclerotic arteries with continuous and pulsed delivery of energy ^[26]. They found that with a prototype ultrasonic wire probe (n = 79 segments), there was gross reduction in vascular lesions as well as microscopic disruption of fibrous and calcified plaques. Normal portions of vessels appeared unaffected by the application of ultrasound. In recent years, experiments also showed the feasibility of inducing HIFU thermal lesions within atheromatous plaques in the femoral arteries of familial hypercholesterolaemic (FH) swine or dog in vivo ^[27,28]. Besides, HIFU begins to be considered as possible methods to treat vulnerable plaque and no significant side effects or adverse events are reported in their studies ^[29,30]. The experimental results indicate that HIFU has significant potential for recanalization of arterial stenosis and complete occlusions.

Current studies of HIFU treatment for arteriosclerosis mainly focus on available experimental models. It should also be appreciated that mathematical models can play a vital role in providing a priori information about the possible outcomes and risks involved to clinical practitioners before the onset of HIFU treatment of atherosclerotic plaque. However, there are few studies on the numerical evaluation of thermal ablation on atherosclerotic plaques. In this paper, a multi-physics field coupling model, which includes HIFU field, thermal field, and blood flow field, is applied to analyze the thermal damage of atherosclerotic plaques. Due to the close proximity of the plaque to the artery, the artery may be suffered from thermal damage by HIFU. To achieve precision in ablating plaque and avoid damage to the artery, accurate positioning of an ultrasound transducer is required. With our method, ultrasound energy produced by a self-focused concave spherical transducer (6.4 cm aperture width, 5.6 cm focal length) is focused at the interface between blood and plaque. The multi-physics systems combing bioheat transfer equation for energy transport equation is applied to predict the evolution of thermal lesion with the Arrhenius equation, which has been widely used in analyzing higher temperature thermal treatments and successfully employed to predict irreversible cell death at hyperthermic temperatures of 43–50 ℃ ^[31]. We focus on Arrhenius scheme for predicting thermal ablation of atherosclerotic plaque on a variety of acoustic factors such as power, frequency, etc.

2.   Mathematical modeling and numerical methods

This study assumes that no cavitation occurs when the maximal pressure does not exceed the threshold pressure in whole computational region. We further hypothesize nonlinear effects and shear waves are neglected when the acoustic wave propagation is linear and also that the amplitude of shear waves in the tissue domain are much smaller than that of the pressure waves.

2.1.   Fluid-solid interactions

We consider fluid-solid interaction coupling an incompressible blood fluid with a hyperelastic solid, i.e., atherosclerotic plaque. The diameters of atherosclerotic plaque and artery lie within the order range of 1–2 cm in the large blood vessels as shown in Figure 1. The description within this section mainly follows ^[32]. By S we denote the Lagrangian reference framework of the solid and by$S\left(t\right)$ its current configuration in Eulerian coordinates. By Ƒ(t) we denote the fluid domain (Including blood and water) at time t matching the solid at the common interface $I\left(t\right)$ = ∂ Ƒ(t)∩$\partial S\left(t\right)$ denote the fluid domain at time $t$ matching the solid at the common interface $\partial I$ = $\partial S\left(t\right)$∩$\partial F\left(t\right)$. By $\mathit \Omega \left(t\right)$ = $\partial {\mathit \Omega }_{S}$∪$\partial {\mathit \Omega }_{F}$∪$\partial I$ we denote the Eulerian fluid-solid interaction domain. The domains $F\left(t\right)$, $S\left(t\right)$ and $\mathit \Omega \left(t\right)$ are all two-dimensional. The boundary of the fluid domain ∂ Ƒ(t) = $I\left(t\right)$∪${\mathit\Gamma }_{f}^{D}\left(t\right)$∪${\mathit\Gamma }_{f}^{out}\left(t\right)$ is split into the interface, a Dirichlet part ${\mathit\Gamma }_{f}^{D}\left(t\right)$ (usually inflow, outflow or rigid walls) and an outflow part ${\mathit\Gamma }_{f}^{out}\left(t\right)$.

For simplicity we assume Dirichlet conditions at the solid boundary apart from the interface $\partial S\left(t\right)$ = $I\left(t\right)$∪${\mathit\Gamma }_{s}^{D}\left(t\right)$. Finally, by I = [0, T] we denote the time interval. With the density ${\rho }_{f }$ the incompressible Navier-Stokes equations are given by

with the right hand side field ${\overrightarrow{f}}_{f}$, boundary data ${\overrightarrow{\upsilon }}_{f}^{D}\left(t\right)$ and the interface velocity ${d}_{t}{\overrightarrow{u}}_{s}$ coming from the coupling to the solid equation. By ${\overrightarrow{\upsilon }}_{f}^{0}$ we denote the initial velocity. The solid problem in term is given on $S$ as

Here we denote by $\overrightarrow{F} = I+\nabla {\overrightarrow{u}}_{S}$ the deformation gradient, by ${\rho }_{S}$ the reference density, by ${\overrightarrow{f}}_{S}$ the right hand side vector, boundary data by ${\overrightarrow{u}}_{S}^{D}$ and by ${\sigma }_{f}({\overrightarrow{\upsilon }}_{f}, {p}_{f}){\overrightarrow{n}}_{f}$ the normal stresses from the coupling to the fluid equations. The attribution of the kinematic interface condition ${\overrightarrow{\upsilon }}_{f}$ = ${d}_{t}{\overrightarrow{u}}_{S}$ to the fluid problem and the dynamic condition to the solid problem is artificial, as the coupled system of (1) and (2) must be considered as one entity. Finally, as material models we consider a Newtonian fluid and a St.Venant Kirchhoff solid

where we denote by $\overrightarrow{E} = \frac{1}{2}({\overrightarrow{F}}_{S}^{T}{\overrightarrow{F}}_{S}-I)$ the Green-Lagrange strain tensor and by ${\nu }_{f}$ the kinematic viscosity. The ${\mu }_{S}$ and ${\lambda }_{S}$ are the Lame′s parameters, Eqs 1, 2 and 3 denote the fluid-solid interaction.

2.2.   Acoustic model

The model uses the pressure acoustics, frequency domain interface to model the acoustic field in the fluid and the tissue domain to obtain the acoustic intensity distribution. The wave equation solved is the homogeneous Helmholtz equation in 2D axisymmetric cylindrical coordinates:

Here $r$ and $z$ are the radial and axial coordinates, $p$ is the acoustic pressure, ${p}_{b}$ is the blood pressure without acoustic field, ${q}_{d}$ is the acoustic dipole domain source (if applicable), $\omega$ is the angular frequency, and $\alpha$ is the frequency dependent in absorptive media. The density, ${ \rho }_{c}$, and the speed of sound, ${c}_{c}$, are complex-valued to account for the material's damping properties.

At the boundary between the solid and fluid domains, denoted ∂Ω_SF, the condition on the boundaries reads

where ${\overrightarrow{u}}_{tt}$ is the solid acceleration, and ${\overrightarrow{F}}_{S}$ is the load (force per unit area) experienced by the solid.

2.3.   Thermal model of radiated tissue

The heat source ${Q}_{h}$ for thermal simulation, given in the plane-wave limit, is then calculated as:

where ${\alpha }_{ABS}$ is the acoustic absorption coefficient, I_h is the acoustic intensity magnitude. Once the pressure is determined by Eqs 3 and 4, the heat transfer in tissue can be obtained by applying Pennes' bioheat transfer equation (BHTE) ^[33]:

where $\rho$, $C$, and $T$ are, respectively, the density, specific heat, and temperature of tissue. Also, ${Q}_{h}$ is the rate of ultrasonic energy deposition per unit volume defined by Eq 6, $w$ is the perfusion rate of blood, ${T}_{0}$ is the initial temperature and equals to 37 ℃. The heat conduction in Eq 7 is based on Fourier's Law of heat conduction and is given as:

where $q$ denotes heat flux, $K$, $\nabla {T}_{t}$, and $\overrightarrow{r}$ are the thermal conductivity, temperature gradient, and position vector, respectively. Considering the non-Fourier behavior of heat conduction, TWMBT can be expressed as follows ^[34,35]

where $\tau$ is thermal relaxation time, which denotes a time lag between heat flux and temperature gradient, leading to significant non-Fourier thermal behavior. Based on (7) and (6), TWMBT can be expressed as below: ^[36]

In the region with a large blood vessel resulting in the local cooling, the energy transport equation is given in the following form ^[33]:

In Eq 11 ${\rho }_{b}{C}_{b}({\overrightarrow{u}}_{b}\cdot \nabla T)$ is related to the convective cooling of vessels and $\overrightarrow{u}$ is the blood velocity. ${\rho }_{b}$, ${C}_{b}$, and ${K}_{b}$ are, respectively, the density, specific heat and thermal conductivity of blood.

2.4.   Thermal damage model

The Arrhenius integral equation for thermal damage is not only related to the characteristics of the biological tissue itself, but also to the temperature T of the tissue at that point and the duration time t at the corresponding temperature. The damage degree of biological tissue is usually estimated by the tissue damage rate D(t) expressed in the form ^[37]

Here, A (1/s) is the material parameter (frequency factor), Ea (J·mol⁻¹) is an activation energy for irreversible damage reaction, T is the thermodynamic temperature, and R ( = 8.314 J·mol⁻¹·K⁻¹) is the universal gas constant. Parameters A and Ea for thermal lesion of atherosclerotic plaques are assumed to be 5.6 × 10⁶³ s⁻¹ and 3.95 × 10⁵ J·mol^{−1 [38]}, respectively. The critical values, D = 0.53, signifies the point when thermal necrosis occurs, and D = 1.0 corresponds to a 63% cell ablation volume ^[39].

2.5.   Numerical simulation

In the current work, the two-dimensional acoustic equation and heat transfer equations along with initial and boundary conditions are solved by utilizing the finite element method (FEM) via COMSOL Multiphysics 5.3 (Pressure Acoustics Frequency Domain, Laminar Flow, Bioheat Transfer, Structural Mechanics; Comsol, Inc.; Burlington, MA). In the model as depicted in Figure 1, unit size free triangle mesh is used to partition the model. A maximum mesh size of λ/6 (λ is the wavelength) in the focal region and a maximum mesh size of λ/4 in the rest of the domain is used for grid generation. The vessel walls are consisted of 5086 vertices of mesh and 9359 triangular elements, and the plaque is built of 547 vertices and 1018 triangular elements. In addition, the blood contains 6944 vertices and 13423 triangular units. In the FEM simulations, the blood–plaque interface and the interface between blood and wall obey boundary conditions given by Eqs 2a, d and are solved by using a partitioned approach, which have been widely applied in solving FSI problems ^{[40,41,42,43]}. In order to estimate the effect of thermal history on plaque lesion evolution, the top and bottom walls of vessel and plaque are modeled as fixed wall boundaries. The Navier-Stokes equations are solved on a freely moving deformed mesh, which constitutes the fluid domain. Finally, the system of governing equations along with initial and boundary conditions are solved using MUMPS solver to obtain the acoustic pressure field. In the analysis of velocity and temperature fields, elements with the Lagrange linear shape functions are used. The system of governing equations along with initial and boundary conditions for the velocity and temperature fields are solved by using MUMPS and GMRES solvers, respectively. A conjugate heat transfer approach is used to analyze thermal interactions between fluids and solids.

3.   Results and discussion

A 1.1-MHz and 15-Wcontinuous ultrasound is used for heating, and other acoustic parameters and thermal parameters are listed in Table 1 ^[44,45,46].

Figure 2 shows the configurations of the predicted pressure and temperature along the axial and lateral directions for a spherically focused transducer with aperture radius of 32 mm at a frequency of 1.1 MHz and for output power ${P}_{out}$ = 15 W at 20 s. The temperature and pressure distributions are strongly dependent on the position relative to the ultrasound focal spot. During the deposition of the focal domain, the peak pressure is 1.48 MPa and the peak temperature is 49.5 ℃. For ensuring continuity in heat flux and temperature at the interface, conjugate heat transfer effects of thermal conduction in the plaque and convection in the blood are considered in the analysis. The temperature elevation at the edge of lower vessel wall near water shows a marked change because of two-way fluid-thermal-solid coupling. It can be observed from the temperature curves in Figure 3 there is no evidence of HIFU-induced temperature damage in the targeted healthy. It seems that the blood plays a role of coolant.

The curves show there is a relative higher temperature in the range from −10 mm to 10 mm along vessel walls. The temperature increases with heating time along the upper vessel wall and the maximum temperature reaches 42.1 ℃ at t = 25 s. The temperature in the lower vessel do not dramatically change in the same range, where the maximum temperature is 38.8 ℃. The peak temperature increases at firstly and then decreases due to the cooling effect of water, which keeps a constant temperature of 25 ℃.

Figure 4 demonstrates the lesion evolution of atherosclerotic plaque with Arrhenius equation. The interval injury integral distributions for selected time steps are presented. The white zone in these figures refers to the values of the lesion integral below 0.5 (undamaged tissue) ^[37], the blue zone refers to the values 0.5 ≤ D < 1.0, it is a partial lesion area, and the red zone (D ≥ 1.0) illustrates the area in which the Arrhenius integral achieved the criterion of tissue ablation. As expected, lesion areas characterize the type of elliptical distribution in the focal region during ultrasound insonation. The thermal lesions are sensitive to variation of temperature and increase with heating time in the domain considered. Figure 5 illustrates the lesion contour for the time of 20 s at different powers. The lesion areas of denaturation are 0 at 10 W, 7.8 mm² at 15 W, and 25.3 mm² at 20 W, respectively, while the ablation areas are 0, 3.6 mm², and 16.7 mm², respectively. Plaque lesion increases obviously with power.

Figure 6 illustrates lesion rate response curves along the axis of the central ultrasound beam for the excitation frequencies of 1.0 MHz, 1.1 MHz and 1.2 MHz. The segment lengths along the center of HIFU beam for 43 mm ≤ r ≤ 48 mm, 48 mm < r ≤ 55 mm, 55 mm < r ≤ 63 mm, and 63 mm < r ≤ 68 mm represent lower vessel wall, blood, plaque, and upper vessel wall, respectively. Not surprisingly, the lesion rate change showed the similar situation. Lesion within whole blood and lower vessel wall cannot be found at 15 W. However, lesion in plaque shows significant difference at different frequencies. The peak lesion rate D increases from 0.75 at 1.0 MHz to 3.8 at 1.1 MHz, and reaches 8.2 at 1.2 MHz. The lesion area of plaque increases obviously with frequency.

In Figure7, we present the thermal lesion region in space𝑟, 𝑧∈[43 68] with different thermal relaxation time. It can be seen the peak lesion rates are 53.9 for 𝜏 = 0, 32.5 for 𝜏 = 10 s, and 17.1 for 𝜏 = 20 s, respectively. At the same time, the peak temperatures obtained with non-Fourier heat conduction equation are 51.4 ℃, 50.6 ℃, and 49.6 ℃, no significant change in the peak temperatures. The results show that the thermal dosage with finite value of relaxation time gives lesser heat-affected area compared to the dosage calculated by considering the infinite propagation of thermal front as in the case of the conventional Fourier model ^[47,48]. For non-homogenous biological samples, non-Fourier heat conduction models work better ^[47,48,49].

4.   Conclusions

In this paper we developed a Fluid-solid thermal model to investigate the response of the thermal lesion of atherosclerotic plaque, subjected to HIFU continuous wave. Our results demonstrate that the model can predict the damage evolution of plaque ablation width and depth. The results also show the model provides comparisons of damage response for different acoustic parameters. It is noted that the results of this study are restricted to the early atherosclerotic plaques, and cavitation effects are not taken into account in this study because the maximum pressure at focus is 1.48 MPa far less than the reported cavitation threshold of 5.1 MPa ^[44]. From the thermal effect point of view, HIFU may present a novel option for ablation of atherosclerotic plaques.

Acknowledgements

Support for this work was supported by the National Nature Science Foundation of China (No.11174077 and No.11474090), Hunan Provincial Nature Science Foundation of China (No. 13JJ3076), the Science Research Program of Education Department of Hunan Province (No. 14A127)

Conflict of interest

The authors declared that they have no conflicts of interest to this work.