Finite Element Simulation of In-Stent Restenosis with Tissue Growth Model ()
1. INTRODUCTION
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 [1]. Coronary artery disease (CAD) is the most common type of heart disease. According to CDC, it claimed 365,914 lives [2] in 2017. Coronary stenting has become a major treatment for coronary stenosis, narrowing a blood vessel leading to restricted blood flow, since its first use in 1986. In the US, 1.2 million patients undergo stent implantations each year [3]. Unfortunately, restenosis, a recurrence of stenosis, often occurs after months or years after the stent implantation. According to [4], nearly one-third of the patients who receive stent implantation require further intervention within six months to reopen stented arteries.
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 [5]. Numerous papers have been published to study the interaction between the stent and the artery. On the one hand, clinical studies are focused on the biological understanding and remedies of restenosis. It has been shown that the biological cause of restenosis could be smooth muscle cell production, endothelial cell abundance, fibroblast growth factors, and so on [6 - 8]. It is observed in in vivo experiments [9] that the vascular wall stress increased immediately after balloon angioplasty. It was also reported that the restenosis rate is closely related to the design of the stents [10]. It is now well-accepted that artery growth is stress-modulated. Therefore, an in-depth study of the stress field in the artery after stent placement becomes necessary. Numerical simulation is undoubtedly a suitable approach to do so.
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 [11]. The stent expansion process was not modeled in this simulation, which failed to capture the contacting shear force that an expanding stent could impose on an artery. A similar study was done in [3] with more realistic settings, and the von Mises stress produced by different stents were compared. Some studies were also more focused on fluids flow impacted by the stents [12 - 14], using fluid-structure interaction methods such as the immersed finite element method [15 - 17].
A novel metal stent is modeled as rigid body in [18] where the simulation is mainly focused on the fluid dynamics in the thoracic aorta. A complete stenting system of balloon-stent-plaque-artery was modeled, and the recoil ratio was carefully examined in [19]. An even more realistic stenting process, including the stent insertion and balloon inflation, was modeled in [20]. Except for the idealized models, patient-specific simulations were also reported. The deployment of Bx Velocity stent in the real human coronary artery was simulated in [21], and the influence of the stent thickness was investigated. The effects of six different stents on a patient-specific artery were compared in [22], showing that the numerical simulation can help select stent selection.
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 [25] further developed the mixture theory by applying it to flow-induced arterial adaptions. Although it has some successful applications, this approach is quite complicated in implementation and lacks generality. A recent two-dimensional continuum mathematical model of restenosis is proposed in [26] in which diffusion-reaction equations are used for modeling the mass balance between species in the arterial wall.
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 [29]. This method has been applied in modeling multiple biological structures including arterial wall [30], heart [31], skin [32], brain [33] etc. For relatively new review in this approach, see [34 , 35].
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.
2. GROWTH MODEL
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 [36].
2.1. Kinematics
The key idea of the kinematics of the growth model is the multiplicative decomposition of deformation gradient
into an elastic part
and a growth part
. Based on the work in [29 , 37], an intermediate configuration
is introduced between the reference configuration
and the current configuration
. As shown in Figure 1,
is the mapping from
to
, and the deformation gradient
is defined in regular way as
. We assume there exists a decomposition:
(1)
Similarly, the right Cauchy-Green tensor
, velocity gradient
, Piola-Kirchhoff stress
and the second Piola-Kirchhoff (PK2) stress
can be decomposed into elastic and growth parts as well. Let
denote the density in the reference configuration,
and
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
,
and
, respectively. In analogy to
, define the Jacobians
and
. The transformations of the volumes become:
(2)
and
(3)
Based on the definitions, we have the following density transformations:
(4)
2.2. Balance Laws
The local form of the balance of mass balances the rate of change of the grown density
with a mass source
:
(5)
Figure 1. The intermediate configuration and the decomposition of deformation gradient.
where the grown density is defined in the reference configuration and can be expressed as
. For soft tissues, it is typically assumed that the density remains constant upon growth, which means
and
. So the mass source can be evaluated as
, where
is the growth part of velocity gradient, and can be easily determined by the growth tensor
.
The growth process is assumed to be quasi-static, therefore the linear momentum balance is simply:
(6)
2.3. Constitutive Equations
An isotropic hyperelastic growth model based on the compressible Neo-Hookean model is used in this study. Denote the Helmholtz free energy as
, the PK2 stress and Kirchhoff stress can be evaluated as:
(7a)
(7b)
where
denotes the elastic part of
. The constitutive moduli
can be obtained by taking the derivative of
with respect to
, and its counterpart in the current configuration can be derived by a push-forward operation:
(8a)
(8b)
The detailed derivation of the hyperelastic stress tensor can be found in [38]. It can be seen immediately from Equations (7) and (8) that both the stress and its moduli differ from their usual form only by an additional variable
, which is determined independently by the governing equation of the growth. Again, this resembles the concept of finite growth resembles finite strain plasticity in the sense that the basic constitutive equations themselves remain unchanged, however, they are evaluated exclusively in terms of the kinematic quantities of the intermediate configuration.
2.4. Growth Function
To complete the stress and constitutive moduli, the growth deformation tensor
must be uniquely determined. It can be either isotropic or anisotropic. In this study, we adopt a simple isotropic form:
(9)
where
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:
(10)
where
is a growth criteria function to enforce a threshold on the growth,
is a weighted function to prevent unlimited growth.
and
are specified with the following forms in [29] with a detailed parameter study:
(11a)
(11b)
where
and
are material constants that are used to control the speed of growth,
is the Mandel stress defined as
, and
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:
(12)
The residual is evaluated as:
(13)
To minimize the residual R, we linearize Equation (13) for tangent moduli K:
(14)
In the local Newton iteration, where the deformation tensor
is known, update the stretch ratio
with
. The corresponding growth part of the deformation tensor
can be calculated using Equation (9). Consequently the elastic part
can be calculated using Equation (1). The elastic stress
and elastic constitutive moduli
can then be obtained using Equations (7) and (8), respectively.
2.5. Baseline Hyperelasticity and Algorithm
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:
(15)
The explicit forms of the stress and constitutive moduli can be obtained based on Equation (15):
(16a)
(16b)
and
(17a)
(17b)
where the operation
and
between two second order tensors
and
are defined as
and
respectively. For more details of the mathematics in con-tinuum mechanics, please refer to [39].
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 [37]:
(18a)
(18b)
and
(19a)
(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
. Algorithm 1 lists the major part of the algorithm.
3. COMPUTATIONAL SETUP
3.1. Geometry
Promus PREMIERTM 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. Figure 2(a) shows the planar sketch of the stent, the dimensions are approximated because the original design is not open for public access. Figure 2(b) shows the planar and 3D design of the stent in millimeter. There are 8 struts in the circumferential direction and 13 layers of struts in the longitudinal direction. Notice that there are 4 connectors between two adjacent layers of struts proximally, while there are only 2 elsewhere. The connectors are arranged helically.
The artery is modeled as a tube with a length of 20 mm, shown in Figure 3. The lumen diameter is 3 mm and the wall thickness is 0.8 mm. The plaque is assumed to have a bell-like shape. It is 9 mm long and
the smallest lumen diameter is 1.8 mm. Therefore, the stenosis rate (measured as
) before
stent operation is 40%.
Figure 2. Geometry of the stent. (a) Planar sketch of the stent; (b) 3D view of the stent.
Figure 3. Geometry of the stenosed artery.
3.2. Boundary Conditions
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 [3], therefore the balloon is neglected in our system and the motion of the stent is simplified as prescribed radial displacement.
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.
Table 1 lists the boundary conditions on the stent and the artery in different steps. Because of the periodic nature of the system, only one half of the stent and artery are modeled and a cylindrical coordinate system is used. It is worth mentioning that the radial displacement in the crimping step is −0.431 mm, which means the stent is crimped onto a cylinder with a diameter of 1.8 mm rather than 1.055 mm. This makes no difference because the stent is in elastic phase during the crimping step, it will result with the same state when it gets in contact with the artery.
,
,
denote the displacements in radial, circumferential, and axial directions respectively.
is measured based on the initial state, not the previous step.
3.3. Materials and Mesh
The major component in the stent is 316LN stainless steel. We model it as an elastoplastic material. The material constants are obtained from [40] as listed in Table 2. The material of the artery is described with the compressible Neo-Hookean as mentioned in Section 2. The material constants and growth parameters are obtained from [30], as listed in Table 3. The material property of the plaque is changing with its development. We assume it is in the beginning of the calcification stage, when it is much stiffer than the healthy artery. Since there is no experimental data for the growth parameters, we adopt them from [29]. This set of parameters has not been calibrated with clinical results; it only describes the growth qualitatively. The growth exponent and adaption speed are relatively high in order to accelerate the growth. The growth criterion is low, which makes growth more evident. Finally the maximum growth factor is set to bound the growth in a reasonable range.
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.
4. RESULTS AND DISCUSSIONS
4.1. Crimping
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. Figure 4 shows the von Mises stress in the artery-stent system at the end of the crimping process. At this point, the artery is stress-free, while the stent is not. From this point on, the stent and the artery are in contact.
4.2. Expanding
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. Figure 5(a) and Figure 5(b) present the geometry of the artery-stent system and the von Mises stress in the artery at the end of the expansion process, respectively.
Table 2. Material constants of the stent.
Table 3. Material constants and growth parameters of the artery.
Figure 4. Von Mises stress in the artery-stent system at the end of crimping process.
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.
4.3. Growth
The growth process starts after the expansion. In our model, the Mandel stress is the source of growth,
Figure 5. Geometry and Von Mises stress at the end of expansion process. (a) Geometry of the artery-stent system; (b) Von Mises stress in the artery.
Figure 6 shows the Mandel stress in the artery during the growth process. As expected, the Mandel stress is gradually decreased because of the growth of the tissue. At the final state, the Mandel stress in the whole artery becomes 0, which indicates the end of the growth. Restenosis is also observed in this figure: the plaque “bounces back” towards the inner direction. The amount of the growth is shown in Figure 7. Starting from the initial state where
, the growth factor increases to the maximum. The growing area expands from the plaque to the vessel wall. The growth in the plaque is much faster than the healthy tissue. However, because the growth is bounded, the lumen diameter of the artery is still larger than the pre-procedure state in Figure 4.
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 Figure 8. At normalized time
, the lumen diameter is pushed to 3 mm. Interestingly, the lumen diameter continues increasing for a short period of time, expanding to 3.2 mm before it starts to decrease. This is due to the compressible nature of the material: in the early stage of the growth, the stress generated from the compression is not large enough to hold the growth. After this point, the restenosis starts to develop, until the growth is terminated by the criterion, which is the so-called biological equilibrium. The final lumen diameter is 2.23 mm, therefore, the loss due to the restenosis is 0.77 mm and the net gain of the stent operation is 0.43 mm. The restenosis rate is 64.2%, which is quite severe but reasonable given the problem setup.
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 Figure 9, the growth is weakened from the inner wall to the outer wall as expected. The relationship between the distance to the inner wall and the growth factor is highly nonlinear. The results show the change in the near-wall area is
Figure 6. Mandel stress in the artery during the growth: (a) The initial state; (b) t = 100; (c) t = 200; (d) t = 300.
Figure 7. Growth factor in the artery during the growth: (a) The initial state; (b) t = 100; (c) t = 200; (d) t = 300.
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.
Figure 9. Growth factors along the radial direction: Point 1 to 5 are distributed uniformly from the inner wall to the outer wall.
5. CONCLUSIONS
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.