^{1}

^{1}

^{*}

In this study, a finite element simulation of in-stent restenosis (ISR) is conducted to simulate the deployment and expansion of a stent in an occluded artery with a contact model and a mechanics-based growth model. A tissue growth model based on the multiplicative decomposition of deformation is applied to investigate the growth of the plaque and artery wall upon the stent’s implantation. Due to the high stresses at the contact points between the stent struts and the tissue, further tissue injury or restenosis is observed. The simulation results show that after the stent deployment, the von Mises stress is significantly larger in the plaque compared to the artery wall, especially in the region that is in contact with the stent. However, the growth of the plaque and artery tends to even out the stress concentration over time. The tissue growth is found to be more significant near the inner wall than the outer layer. A 0.77 mm restenosis is predicted, which agrees with published clinical observations. The features of the artery growth are carefully analyzed, and the underlying mechanism is discussed. This study is the first attempt to apply finite element analysis to artery restenosis, which establishes a framework for predicting ISR’s occurrence and severity. The results also provide insights into understanding the underlying mechanism of in-stent restenosis.

Cardiovascular disease (CVD) is the leading cause of mortality in the United States in men and women of every major ethnic group. It accounted for nearly 840,768 deaths in 2016 [

Restenosis is characterized by neo-intima formation, resulting from proliferating medial vascular smooth muscle cells migrating to the intima, and the synthesis of matrix components, mainly collagen, with resultant stenosis and constrictive tissue remodeling [

In one of the early studies, finite element analysis was carried out to examine the maximum principal stress in the stenosed artery caused by two different types of stents [

A novel metal stent is modeled as rigid body in [

All of the previous simulation studies stopped at the calculation of the stress field in the stented artery. In this study, we attempt to simulate the real growth of the artery under the impact of the stent’s deployment. There are various ways to model mechanical-induced growth; for example, the multi-constituent open systems based on the theory of mixtures [23 , 24]. Humphrey and Rajagopal [

Unlike the mathematical models cited previously, a mechanical growth model based on the multiplicative decomposition of the deformation was developed [27 , 28]. In this open system thermodynamics framework, an incompatible “growth” configuration was postulated between the material configuration and current configuration, resembling the idea of multiplicative decomposition in the context of finite plasticity, see [

Due to its realisms in presenting tissue growth and simplicity in implementation, the growth model based on the incompatible growth configuration and the decomposition of the deformation gradient into an elastic and growth part is adopted in our study to characterize the long term growth of the arterial wall after implanting the stent. The rest of the manuscript is organized as following: in Section 2, we briefly review the mechanical growth model; in Section 3, the model setup, as well as the geometry, mesh, and boundary conditions are introduced; the results are presented in Section 4, followed by discussions; and finally conclusions are drawn in Section 5.

In this section, we briefly review the mechanical growth model developed in [27 , 28]. The model is described in the aspects of kinematics, equilibrium equations and constitutive equations. This growth model was successfully applied to the simulation of vessel tissue remodeling with residual stress [

The key idea of the kinematics of the growth model is the multiplicative decomposition of deformation gradient F into an elastic part F e and a growth part F g . Based on the work in [29 , 37], an intermediate configuration B g is introduced between the reference configuration B 0 and the current configuration B t . As shown in

F = F e ⋅ F g (1)

Similarly, the right Cauchy-Green tensor C , velocity gradient l , Piola-Kirchhoff stress P and the second Piola-Kirchhoff (PK2) stress S can be decomposed into elastic and growth parts as well. Let ρ 0 denote the density in the reference configuration, ρ t and ρ g denote its counterparts in the spatial configuration and the intermediate configuration, respectively. Similarly, define the volume of a material particle in the reference, spatial and intermediate configurations as d V , d v and d V g , respectively. In analogy to J = det ( F ) , define the Jacobians J e = det ( F e ) and J g = det ( F g ) . The transformations of the volumes become:

d v = J d V , d V g = J g d V , d v = J e d V g (2)

and

J = J e J g (3)

Based on the definitions, we have the following density transformations:

ρ g d V g = ρ t d v , ρ g = J e ρ t (4)

The local form of the balance of mass balances the rate of change of the grown density ρ ¯ 0 with a mass source R 0 :

ρ ¯ ˙ 0 = R 0 (5)

where the grown density is defined in the reference configuration and can be expressed as ρ ¯ 0 = J ρ t = J g ρ g . For soft tissues, it is typically assumed that the density remains constant upon growth, which means ρ g = const and ρ ˙ g = 0 . So the mass source can be evaluated as R 0 = ρ ¯ 0 tr L g , where L g is the growth part of velocity gradient, and can be easily determined by the growth tensor F g .

The growth process is assumed to be quasi-static, therefore the linear momentum balance is simply:

div σ = 0 (6)

An isotropic hyperelastic growth model based on the compressible Neo-Hookean model is used in this study. Denote the Helmholtz free energy as ϕ ( C , F g ) , the PK2 stress and Kirchhoff stress can be evaluated as:

S = 2 ∂ Ψ ∂ C = 2 ∂ Ψ ∂ C e : ∂ C e ∂ C = F g − 1 ⋅ S e ⋅ F g − T (7a)

τ = F ⋅ S ⋅ F T = F e ⋅ S e ⋅ F e T (7b)

where S e denotes the elastic part of S . The constitutive moduli ℂ can be obtained by taking the derivative of S with respect to C , and its counterpart in the current configuration can be derived by a push-forward operation:

ℂ = 2 d S ( F , F g ) d C = 2 ∂ S ( F , F g ) ∂ C | F g + 2 [ ∂ S ∂ F g : ∂ F g ∂ θ g ] ⊗ ∂ θ g ∂ C | F (8a)

ℂ = ( F ⊗ ¯ F ) : ℂ : ( F T ⊗ ¯ F T ) (8b)

The detailed derivation of the hyperelastic stress tensor can be found in [

To complete the stress and constitutive moduli, the growth deformation tensor F g must be uniquely determined. It can be either isotropic or anisotropic. In this study, we adopt a simple isotropic form:

F g = θ g I (9)

where θ g denotes an internal variable which is often referred to as growth factor. The growth multiplier is 1 in the plain elastic deformation, greater than 1 in growth, and smaller than 1 in atrophy.

To model the evolution of the growth multipliers, the rate of growth multipliers is usually specified as:

θ ˙ g = k g ( θ g ) ⋅ ϕ g ( F e ) (10)

where ϕ g is a growth criteria function to enforce a threshold on the growth, k g is a weighted function to prevent unlimited growth. ϕ g and k g are specified with the following forms in [

ϕ g = tr ( M e ) − M e crit (11a)

k g = 1 τ ( θ max − θ g θ max − 1 ) γ (11b)

where τ and γ are material constants that are used to control the speed of growth, M e is the Mandel stress defined as M e = C e ⋅ S e , and M e crit is the criteria of growth. A local Newton iteration can be used to solve for the growth factor. Applying backward Euler scheme to Equation (10), we have:

θ ˙ g = ( θ g − θ g n ) / Δ t (12)

The residual is evaluated as:

R = θ g − θ g n − 1 τ ( θ max − θ g θ max − 1 ) γ ( tr ( M e ) − M e crit ) Δ t (13)

To minimize the residual R, we linearize Equation (13) for tangent moduli K:

K = ∂ R ∂ θ g = 1 − ( k g ∂ ϕ g ∂ θ g + ϕ g ∂ k g ∂ θ g ) Δ t (14)

In the local Newton iteration, where the deformation tensor F is known, update the stretch ratio θ g with θ g ← θ g − R / K . The corresponding growth part of the deformation tensor F g can be calculated using Equation (9). Consequently the elastic part F e can be calculated using Equation (1). The elastic stress S e and elastic constitutive moduli ℂ e can then be obtained using Equations (7) and (8), respectively.

In this preliminary work, we use compressible Neo-Hookean as the baseline hyperelastic model since it has been validated in multiple applications [28 - 30 , 37]. But there is no fundamental difficulty preventing other models being used. The compressible Neo-Hookean free energy is written as:

Ψ = λ 2 ln 2 J e + 1 2 μ [ C e : I − 3 − 2 ln J e ] (15)

The explicit forms of the stress and constitutive moduli can be obtained based on Equation (15):

S e = [ λ ln J e − μ ] C e − 1 + μ I (16a)

τ e = ( λ ln J e − μ ) I + μ B e (16b)

and

ℂ e = ( μ − λ ln J e ) ( C e − 1 ⊗ ¯ C e − 1 + C e − 1 ⊗ _ C e − 1 ) + λ C e − 1 ⊗ C e − 1 (17a)

ℂ e = ( μ − λ ln J e ) ( δ i k δ j l + δ i l δ j k ) + λ δ i j δ k l (17b)

where the operation ⊗ ¯ and ⊗ _ between two second order tensors A and B are defined as

( A ⊗ ¯ B ) i j k l = A i k B j l and ( A ⊗ _ B ) i j k l = A i l B j k respectively. For more details of the mathematics in con-tinuum mechanics, please refer to [

For completeness, we list the specific expressions of Equations (7) and (8) using compressible Neo-Hookean as the baseline hyperelasticity, the detailed derivation can be found in [

S = 1 θ e 2 S e (18a)

τ = τ e (18b)

and

ℂ = 1 θ g 4 ℂ e − 4 θ g 5 k g K Δ t ( S e + 1 2 ℂ e : C e ) ⊗ ( 1 2 C e : ℂ e + S e ) (19a)

ℂ = ℂ e + 4 θ g k g K Δ t ( F e ⊗ ¯ F e ) : [ ( S e + 1 2 ℂ e : C e ) ⊗ ( 1 2 C e : ℂ e + S e ) ] : ( F e T ⊗ ¯ F e T ) (19b)

Growth model can be easily implemented within the existing nonlinear finite element framework. In this research the updated Lagrangian formulation is used. Except for the modifications to the Kirchhoff stress and the corresponding constitutive moduli, all the additional computational work it introduces is the iterative solution of the stretch ratio θ g . Algorithm 1 lists the major part of the algorithm.

Promus PREMIER^{TM} is a relatively new coronary stent designed by Boston Scientific that has a higher strength and less restenosis rate compared to older generations. Here, a 16 mm length model is created with SOLIDWORKS (Dassault Systémes). It is constructed from a cylinder base with a diameter of 2.5 cm and crimped to a catheter with a diameter of 1.055 cm. The thickness of the strut is 0.081 mm.

The artery is modeled as a tube with a length of 20 mm, shown in

the smallest lumen diameter is 1.8 mm. Therefore, the stenosis rate (measured as 1 − D blocked lumen D original lumen ) before

stent operation is 40%.

In reality, the stent has to be crimped onto a catheter before inserted into the vessel. After the stent is in place, a folded balloon will be inflated to expand the stent. In this study, we are interested in the interaction between the stent and the artery, besides the expansion of the effective part of the stent is rather uniform [

The simulation is divided into 3 steps: stent crimping, stent expanding, and vessel growth. In the initial stress-free state, the diameter of the stent is larger than the lumen diameter of the stenosed artery. To simulate the crimping process, in the first step, the stent is compressed in the radial direction. The second step is to simulate the expansion of the stent and its interaction with the artery. The stent will expand freely before it reaches the inner wall of the artery. Then the stent will push inner wall of the plaque, expanding the lumen diameter to 3 mm; the last step is the main focus of this study, where we predict the growth of the plaque and artery using the growth model introduced in Section 2.

Stent | Artery | |
---|---|---|

Crimping | U r = − 0.431 mm | - |

U θ = 0 on side planes | ||

U z = 0 on reference point | ||

Expanding | U r = 0.169 mm | U θ = 0 |

U θ = 0 on side planes | ||

U z = 0 on reference point | U z = 0 on top and bottom planes | |

Growth | - | U r = 0 on the outer wall |

U θ = 0 | ||

U z = 0 on top and bottom planes |

U r , U θ , U z denote the displacements in radial, circumferential, and axial directions respectively. U r is measured based on the initial state, not the previous step.

The major component in the stent is 316LN stainless steel. We model it as an elastoplastic material. The material constants are obtained from [

In this analysis, the displacement of the stent is prescribed, therefore it is modeled with shell element to simplify the computation. A total of 7261 elements are used. The artery is modeled with C3D8 element since the geometry is simple and material is compressible, 25,560 elements are used.

The diameter of the outer surface of the stent is 2.662 mm in stress-free state. It should be crimped onto a 1.055 mm-diameter catheter before being inserted into the artery. However, since the crimping process is in the elastic range, we only need to crimp it to 1.8 mm, which is as wide as the stenosed artery, so that it does not penetrate the artery before the expansion begins.

As the stent continues expanding, the plaque will be pushed to the artery wall. The goal of the operation is to recover the lumen diameter to 3 mm.

Young’s modulus (MPa) | 1.93×10^{5} | |
---|---|---|

Poisson’s ratio | 0.33 | |

Plastic behavior | Yield stress (MPa) | Plastic strain |

304 | 0 | |

484 | 0.0928 | |

612 | 0.179 | |

868 | 0.4 |

artery | plaque | |
---|---|---|

μ (MPa) | 1.15385 | 11.5385 |

λ (MPa) | 1.7307 | 17.307 |

growth criterion (MPa) | 0 | 0 |

maximum growth factor | 1.3 | 1.3 |

growth exponent γ | 5 | 5 |

adaption speed 1 / τ | 0.5 | 0.5 |

It can be seen that after the stent operation, the plaque is significantly compressed, and the inner wall of the vessel becomes flat. The von Mises stress in the healthy tissue is negligible compared to the plaque. The contact area can be clearly observed based on the stress field. Maximum stress in the plaque indicates the danger of the rupture.

The growth process starts after the expansion. In our model, the Mandel stress is the source of growth,

To visualize the restenosis quantitively, we look at the change of the lumen diameter from the pre-procedure state to the final stable state after restenosis as shown in

To compare the growth at different locations in the artery, the growth factors at 5 points along the radial direction in the middle of the artery wall are sampled. It can be seen from

not as significant as in the middle. Moreover, the change of the growth factor over the normalized time is also highly nonlinear. This is likely due to the sensitivity of the growth model. Unfortunately only a limited number of parameter studies are done in literature to perform a quantitate analysis. The qualitative findings are found to be consistent with clinical observations.

Restenosis is a common side effect of the stent procedure which is caused by the tissue growth. Clinical study has been carried out for years but the mechanism of the growth is still not clear. Researchers in the field of computational mechanics have only been involved in the area of biological growth recently yet there are already plenty of outcomes. The stent procedure has been studied in many previous papers. However, in this study, we push the study forward by applying a stress-induced growth model to the post-procedure artery. The growth of the artery tissue is investigated. Qualitatively, it is found that the maximum growth factor is mostly concentrated in the area that is in contact with the stent, since it is where the stress is the highest. This implies that the goal of the design of the future stents should be focused on reducing stress concentration. Quantitively, the lumen diameter of the post-procedure artery is found to be re-narrowed from 3 mm to 2.23 mm, which is reasonable based on clinical observations.

It can be seen from this study that finite element analysis is a promising tool to study some of the complicated biological growth that in vivo investigations are difficult to apply to. However, there are also unsolved problems: first, in order to characterize the growth more realistically, experiments must be done to find the proper growth parameters; second, the growth is assumed to be isotropic, which is not necessarily true because we know the vessel is highly anisotropic. More sophisticated growth model may be required to capture the anisotropy.

The authors declare no conflicts of interest regarding the publication of this paper.