This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

11institutetext: ∗ corresponding author

Manjunatha, K., Reese, S. 22institutetext: RWTH Aachen University, Institute of Applied Mechanics, Aachen, Germany
22email: [email protected], [email protected]
33institutetext: Behr, M. 44institutetext: RWTH Aachen University, Chair for Computational Analysis of Technical Systems, Aachen,
Germany
44email: [email protected] 55institutetext: Vogt, F. 66institutetext: RWTH Aachen University, Department of Cardiology, Pulmonology, Intensive Care and Vascular Medicine, Aachen, Germany
66email: [email protected]

Submitted for publication in a “Festschrift” for Peter Wriggers on December 1,2020.

Finite element modelling of in-stent restenosis

Kiran Manjunatha∗    Marek Behr    Felix Vogt    Stefanie Reese
Abstract

From the perspective of coronary heart disease, the development of stents has come significantly far in reducing the associated mortality rate, drug-eluting stents being the epitome of innovative and effective solutions. Within this work, the intricate process of in-stent restenosis is modelled considering one of the significant growth factors and its effect on constituents of the arterial wall. A multiphysical modelling approach is adopted in this regard. Experimental investigations from the literature have been used to hypothesize the governing equations and the corresponding parameters. A staggered solution strategy is utilised to capture the transport phenomena as well as the growth and remodeling that follows stent implantation. The model herein developed serves as a tool to predict in-stent restenosis depending on the endothelial injury sustained and the protuberance of stents into the lumen of the arteries.

Keywords: in-stent restenosis, smooth mucsle cells, platelet-derived growth factor, extracellular matrix, growth

1 Introduction

Coronary heart disease (CHD) is one of the most prominent causes of mortality and affected 126.5 million lives worldwide as of 2017. CHD is characterised by constriction of the coronary artery arising due to the build-up of plaque which consists of lipids, calcifications, and macrophages. This condition, termed atherosclerosis, leads to the narrowing of pathways for blood flow as well as the loss of elasticity of the arterial wall. Percutaneous coronary intervention (PCI) is the process of placing reinforcement structures called stents within these constricted blood vessels to normalize the blood flow. Unfortunately, this interventional procedure is associated with the risk of in-stent restenosis and stent thrombosis. Drug-eluting-stents have emerged in recent years as the most viable option for revascularisation, significantly reducing the risks associated with PCI. But the mechanisms involved in the restenotic process remain incompletely understood. Quantifying the probability and level of restenosis via simulation of the underlying mechanisms will help in addressing the risks more precisely.

Several computational approaches have been developed in this regard. A multi-scale framework involving finite element models and agent-based models has been conceptualised in zahedmanesh2012 , where the restenotic process was unidirectionally coupled to stresses on the arterial walls. A more recent development on similar grounds has been made in li2019 , incorporating bidirectional coupling between finite element and agent-based models. A general continuum-based model describing growth and remodelling of soft tissues humphrey2002 , considering the evolution of constituents in the arterial wall using the constrained-mixture theory, also serves as a computational tool for restenosis prediction. With temporal averaging, a homogenised constrained-mixture model has been developed in cyron2016 as an extension to the constrained-mixture theory. Alternatively, highly resolved transport phenomena occurring in the arterial wall have also been utilised to model the pathophysiology of vascular diseases. Venous neointimal hyperplasia has been quantified considering the inherent pathogenic mechanisms involving growth factors in budu2008 . A coupled multi-physical approach for quantifying atherosclerosis has been described in thon2018 , wherein the mechanics of bloodflow, mass transport and arterial wall mechanics have been unified under a single framework. An in-stent restenosis predictive model developed on a similar foundation incorporating damage on the arterial wall can be found in escuer2019 . An analogous attempt is made hereof for the prediction and quantification of the after-effects of stent implantation.

1.1 Structure of the arterial wall

The arterial wall is composed of three layers: intima, media, and adventitia. The region where the blood flow occurs is called the lumen. Intima is the innermost layer of the arterial wall, immediately adjacent to the lumen. It is consisted of a monolayer of endothelial cells accompanied by a few layers of smooth muscle cells and extends up to the internal elastic lamina. The media occupies the space between the internal and external laminae and is mainly composed of collagen, elastin and smooth muscle cells. Adventitia, the outermost layer, contains loose connective tissue composed of elastin, collagen and fibroblasts.

1.2 In-stent restenosis and platelet-derived growth factor

In-stent restenosis refers to the accumulation of new tissue within the intima leading to a diminished cross-section of the lumen post stent implantation. The underlying mechanism is called neointimal hyperplasia. It is a collaborative effect of migration and proliferation of smooth muscle cells (SMC) in the arterial wall, regulated by intricate signalling pathways that are triggered by certain stimuli, either internal or external to the arterial wall.

The platelet-derived growth factor (PDGF) is a dimer of two peptides linked by a disulfide bond. It has been implicated in vascular remodelling processes, including neointimal hyperplasia, that follow an injury to arterial wall koyama1994 . This can be attributed to its mitogenic and chemoattractant properties. PDGF is secreted by an array of cellular species namely the endothelial cells, SMC, fibroblasts, macrophages and platelets.

The stent implantation procedure damages the endothelial monolayer. Also, depending on the arterial overstretch achieved during the implantation, stent struts partially obstruct the blood flow creating vortices in their wake regions. This causes oscillatory wall shear stresses and hence further damages to the endothelium koskinas2012 . PDGF, which is stored in the alpha-granules of the aggregated platelets at endothelial injury sites, is released into the arterial wall. The presence of PDGF upregulates matrix metalloproteinases (MMP) production in the arterial wall. Extracellular matrix (ECM) is a network of collagen and glycoproteins surrounding the SMC and are degraded in the presence of MMP. SMC, which are usually held stationary by the ECM, are rendered free for migration under the action of MMP. Also, a degraded ECM encourages the proliferation of SMC under the presence of growth factors. The focal adhesion sites created due to cleaved proteins in the ECM assist in the migration of SMC. The direction of migration is influenced by the number of adjacent focal adhesion sites available and the local concentration gradient in PDGF. This directional movement of SMC is termed chemotaxis, and results in the accumulation of the proliferated and migrated SMC in the intima of the arterial wall. A positive feedback loop might occur wherein the migrated SMC create a further obstruction in the blood flow and subsequent upregulation of PDGF. The uncontrolled growth of vascular tissue that follows will eventually lead to a complete blockage of the lumen.

2 Mathematical modelling

In-stent restenosis is understood to be a convoluted phenomenon involving multiple constituents. Only the behaviours of key constituents identified from the perspective of the processes described in Sect.1.2 are modelled. Blood flow in the lumen is not included in the model and the arterial wall is considered to be made of only the intima and media layers.

2.1 Transport phenomena

The transport of constituents within the arterial wall is governed by a set of advection-reaction-diffusion equations, the general structure of which for a scalar field ϕ\phi is

∂ϕ∂t=−∇⋅(ϕ​𝒗)⏟advection+∇⋅(k​∇ϕ)⏟diffusion+Ts⏟source−Tr⏟reaction,\displaystyle{\frac{\partial\phi}{\partial t}}=-\underbrace{\nabla\cdot\left(\phi\,\boldsymbol{v}\right)}_{\text{advection}}+\underbrace{\nabla\cdot\left(k\,\nabla\phi\right)}_{\text{diffusion}}+\underbrace{T_{s}}_{\text{source}}-\underbrace{T_{r}}_{\text{reaction}}, (1)

where 𝒗\boldsymbol{v} denotes the velocity of the medium of transport and kk, the diffusivity of ϕ\phi in the medium. To mathematically model in-stent restenosis, PDGF, ECM and SMC are considered to be the key ingredients. Hence the general structure above is adapted to reflect their respective behaviours. The arterial wall is considered to be quasi-static and hence advective terms arising from the movement of the wall are ignored (𝒗≈𝟎)(\boldsymbol{v}\approx\boldsymbol{0}). The symbols and units associated with the transport quantities are declared here for clarity.

Table 1: Transport variables
Variable Symbol Units
\svhline PDGF concentration cPc_{{}_{P}} mol/mm3
ECM density ρE\rho_{{}_{E}} mol/mm3
SMC density ρS\rho_{{}_{S}} cells/mm3

Platelet-derived growth factor

PDGF is assumed to diffuse throughout the arterial wall. Also, it is internalised by the SMCs during their migration and proliferation. The evolution of the concentration of PDGF is thus modelled by the following equation:

∂cP∂t=∇⋅(DP​∇cP)⏟diffusion​−α​ρS​cP⏟consumption by SMC,\displaystyle{\frac{\partial c_{{}_{P}}}{\partial t}}=\underbrace{\nabla\cdot\left(D_{{}_{P}}\,\nabla c_{{}_{P}}\right)}_{\text{diffusion}}\underbrace{-\,\,\alpha\,\rho_{{}_{S}}\,c_{{}_{P}}}_{\text{consumption by SMC}}, (2)

where α\alpha refers to an internalisation coefficient and ρS\rho_{{}_{S}}, the density of SMCs. In comparison to Eq. 1, we see that only the terms associated with the diffusive and the reactive phenomena appear here. Although PDGF is produced by other cellular species, the net effect observed is degradation. The diffusion coefficient DPD_{{}_{P}} can theoretically vary between positions in the arterial wall but is considered to be constant here.

Extracellular matrix

The major constituent of ECM is collagen, which is considered to be non-diffusive. SMCs recognize a degraded ECM and synthesize collagen. A source term, in the form of a logistic function, is introduced in this regard and an asymptotic threshold for collagen density ρE,t​h\rho_{{}_{E,th}} prescribed. Collagen is degraded by MMP, which in turn is regulated by PDGF. A reactive term is introduced to take care of this degradation. The evolution of ECM density hence reads as follows:

∂ρE∂t=β​ρS​(1−ρEρE,t​h)⏟synthesis by SMC​−γ​cP​ρE⏟degradation due to MMP.\displaystyle{\frac{\partial\rho_{{}_{E}}}{\partial t}}=\underbrace{\beta\,\rho_{{}_{S}}\,\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)}_{\text{synthesis by SMC}}\underbrace{-\,\,\gamma\,c_{{}_{P}}\,\rho_{{}_{E}}}_{\text{degradation due to MMP}}. (3)

Here, β\beta and γ\gamma refer to synthesis and degradation rate coefficients respectively.

Smooth muscle cells

The migration of SMC under the influence of PDGF is modelled using a chemotaxis term keller1971 , which can be interpreted as a pseudo-advective term wherein the velocity of the medium is replaced by the gradient in the ECM density. This corresponds to the fact that higher degradation in the ECM results in higher focal adhesion sites for SMC migration. Additionally, a source term is introduced to model the proliferation of SMC under the presence of PDGF. Hence the evolution of the SMC density is prescribed using the equation:

∂ρS∂t=∇⋅(χ​cP​(1−ρEρE,t​h)​ρS​∇ρE)⏟chemotaxis due to ECM degradation+κ​cP​(1−ρEρE,t​h)​ρS⏟proliferation due to PDGF.\displaystyle{\frac{\partial\rho_{{}_{S}}}{\partial t}}=\underbrace{\nabla\cdot\left(\chi\,c_{{}_{P}}\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)\,\rho_{{}_{S}}\,\nabla\rho_{{}_{E}}\right)}_{\text{chemotaxis due to ECM degradation}}+\underbrace{\kappa\,c_{{}_{P}}\,\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)\,\rho_{{}_{S}}}_{\text{proliferation due to PDGF}}. (4)

The chemotactic sensitivity χ\chi and the proliferation constant κ\kappa are scaled by the concentration of PDGF as well as the logistic coefficient dependent on the ECM density. This takes care of the fact that the migration and proliferation effects increase with an increase in ECM degradation, in the presence of PDGF.

2.2 Arterial wall mechanics

The continuum mechanical description of in-stent restenosis requires first the definition of the kinematics of growth and remodelling of the arterial wall, incorporating relevant quantities from the transport phenomena previously described. Constitutive equations hence based on the kinematics, in combination with the balance of linear momentum, will complete the mathematical model description.

2.3 Kinematics

A deformation map 𝝋\boldsymbol{\varphi} of a particle at position 𝑿\boldsymbol{X} in the reference configuration Ω0\Omega_{0} at time t0t_{0} to its position 𝒙\boldsymbol{x} in the current configuration Ω\Omega at time tt is provided by the deformation gradient 𝑭=∇𝑿𝝋​(𝑿,t)\boldsymbol{F}=\nabla_{{}_{\boldsymbol{X}}}\boldsymbol{\varphi}(\boldsymbol{X},t). The right Cauchy-Green tensor is further defined as 𝑪=𝑭T​𝑭\boldsymbol{C}=\boldsymbol{F}^{T}\,\boldsymbol{F}.

For a description of growth, a multiplicative decomposition of the deformation gradient is adopted. Assumption of an intermediate incompatible configuration which achieves a locally stress-free state forms the basis of this split himpel2005 ; rodriguez1994 . An additional elastic deformation is needed to ensure the compatibility of the total deformation. Hence the total deformation gradient will be 𝑭=𝑭e​𝑭g\boldsymbol{F}=\boldsymbol{F}_{e}\,\boldsymbol{F}_{g} and the elastic right Cauchy-Green tensor reads 𝑪e=𝑭eT​𝑭e\boldsymbol{C}_{e}=\boldsymbol{F}_{e}^{T}\,\boldsymbol{F}_{e}.

The growth deformation gradient is specified as 𝑭g=ϑ​𝑰\boldsymbol{F}_{g}=\vartheta\,\boldsymbol{I}, under the assumption of isotropic growth. Hence the compatible elastic deformation gradient will be 𝑭e=ϑ−1​𝑭.\boldsymbol{F}_{e}=\vartheta^{-1}\,\boldsymbol{F}.

To calculate the growth stretch ϑ\vartheta, it is hypothesised that an increase in volume at a point due to growth at any time tt is proportional to the additional SMC transported to that point. A similar idea is postulated for aggregation of foam cells in thon2018 . If vsv_{s} is the volume occupied by a single SMC and NsN_{s} the number of additional SMC at a point, then Δ​Vg=vs​Ns\Delta V_{g}=v_{s}\,N_{s}. vsv_{s} can be prescribed via the definition of the SMC density of a healthy artery ρS,h\rho_{{}_{S,h}}, using vs=(ρS,h)−1v_{s}=(\rho_{{}_{S,h}})^{-1}. If VgV_{g} is the volume in the intermediate growth configuration Ωg\Omega_{g} at time tt and V0V_{0} is that in the reference configuration Ω0\Omega_{0}, then we can arrive at an expression for ϑ\vartheta through the following exercise.

Δ​Vg=Vg−V0\displaystyle\Delta V_{g}=V_{g}-V_{0} =\displaystyle= (ρS,h)−1​Ns\displaystyle(\rho_{{}_{S,h}})^{-1}\,N_{s} (5)
∫Ωg1​𝑑vg−∫Ω01​𝑑V\displaystyle\displaystyle{\int_{\Omega_{g}}}1\,dv_{g}-\displaystyle{\int_{\Omega_{0}}}1\,dV =\displaystyle= (ρS,h)−1​∫Ω(ρS−ρS,h)​𝑑v\displaystyle(\rho_{{}_{S,h}})^{-1}\,\displaystyle{\int_{\Omega}}(\rho_{{}_{S}}-\rho_{{}_{S,h}})\,dv (6)

Pulling back the quantities to the initial configuration,

∫Ω0Jg​𝑑V−∫Ω01​𝑑V=(ρS,h)−1​∫Ω0J​(ρS−ρS,h)​𝑑V\displaystyle{\int_{\Omega_{0}}}J_{g}\,dV-\displaystyle{\int_{\Omega_{0}}}1\,dV=(\rho_{{}_{S,h}})^{-1}\,\displaystyle{\int_{\Omega_{0}}}J\,(\rho_{{}_{S}}-\rho_{{}_{S,h}})\,dV (7)

Since Eq. 7 has to hold true locally,

Jg−1\displaystyle J_{g}-1 =\displaystyle= (ρS,h)−1​J​(ρS−ρS,h)\displaystyle(\rho_{{}_{S,h}})^{-1}\,J\,(\rho_{{}_{S}}-\rho_{{}_{S,h}}) (8)
⟹Jg\displaystyle\implies J_{g} =\displaystyle= 1+J​(ρSρS,h−1).\displaystyle 1+\displaystyle{J\,\left(\frac{\rho_{{}_{S}}}{\rho_{{}_{S,h}}}-1\right)}. (9)

If dd is the dimensionality of the problem, then Jg=det(𝑭g)=ϑdJ_{g}=\det(\boldsymbol{F}_{g})=\vartheta^{d}. Therefore,

ϑ=(1+J​(ρSρS,h−1))1/d.\vartheta=\left(1+\displaystyle{J\,\left(\frac{\rho_{{}_{S}}}{\rho_{{}_{S,h}}}-1\right)}\right)^{1/d}. (10)

2.4 Hyperelastic constitutive model

The constitutive model is derived from the hyperelastic models presented in gasser2006 ; holzapfel2000 ; nolan2014 , wherein the arterial wall is assumed to be composed of two helices of collagen fibres, with respective helix angles θia\theta^{a}_{i}, embedded in an isotropic ground substance. The associated Helmholtz free energy per unit volume is hence split into an isotropic and an anisotropic part, both of which are dependent on the elastic right Cauchy-Green tensor.

ψ​(𝑪e,𝑯1,𝑯2)\displaystyle\psi(\boldsymbol{C}_{e},\boldsymbol{H}_{1},\boldsymbol{H}_{2}) =\displaystyle= ψi​s​o​(𝑪e)+ψa​n​i​(𝑪e,𝑯1,𝑯2)\displaystyle\psi_{iso}(\boldsymbol{C}_{e})+\psi_{ani}(\boldsymbol{C}_{e},\boldsymbol{H}_{1},\boldsymbol{H}_{2}) (11)
ψi​s​o​(𝑪e)\displaystyle\psi_{iso}(\boldsymbol{C}_{e}) =\displaystyle= μa2​(tr​𝑪e−3)−μa​ln​Je+Λa4​(Je2−1−2​ln​Je)\displaystyle\displaystyle{\frac{\mu^{a}}{2}}\left(\text{tr}\,\boldsymbol{C}_{e}-3\right)-\mu^{a}\,\text{ln}\,J_{e}+\displaystyle{\frac{\Lambda^{a}}{4}}\left(J_{e}^{2}-1-2\,\text{ln}\,J_{e}\right) (12)
ψa​n​i​(𝑪e,𝑯1,𝑯2)\displaystyle\psi_{ani}(\boldsymbol{C}_{e},\boldsymbol{H}_{1},\boldsymbol{H}_{2}) =\displaystyle= k12​k2​∑i=1,2(exp​[k2​⟨Ei⟩2]−1)\displaystyle\displaystyle{\frac{k_{1}}{2k_{2}}}\sum_{i=1,2}\left(\text{exp}\left[k_{2}\langle E_{i}\rangle^{2}\right]-1\right) (13)

This stems from the hypothesis that in the stress-free incompatible growth configuration, even the collagen fibres grow due to the synthesis by SMC and hence no residual stresses arise in the incompatible state. The dependence of ψ\psi on 𝑪e\boldsymbol{C}_{e} embeds its dependence on 𝑭g\boldsymbol{F}_{g} via

𝑪e=𝑭g−T​𝑪​𝑭g−1=ϑ−2​𝑪.\boldsymbol{C}_{e}=\boldsymbol{F}_{g}^{-T}\,\boldsymbol{C}\,\boldsymbol{F}_{g}^{-1}=\vartheta^{-2}\boldsymbol{C}. (14)

The generalised structure tensors 𝑯i\boldsymbol{H}_{i} are constructed based on the local collagen orientations in the reference configuration 𝒂0​i\boldsymbol{a}_{0i} using

𝑯i=κa​𝑰+(1−3​κa)​𝒂0​i⊗𝒂0​i,\boldsymbol{H}_{i}=\kappa^{a}\,\boldsymbol{I}+\left(1-3\,\kappa^{a}\right)\,\boldsymbol{a}_{0i}\otimes\boldsymbol{a}_{0i}\,\,, (15)

where κa\kappa^{a} is a dispersion parameter. The Green-Lagrange strain like parameter EiE_{i} is calculated utilising the relationship Ei=𝑯i:𝑪e−1,E_{i}=\boldsymbol{H}_{i}:\boldsymbol{C}_{e}-1, where the definition of scalar product of second order tensors 𝑨:𝑩=Ai​j​Bi​j\boldsymbol{A}:\boldsymbol{B}=A_{ij}\,B_{ij} (Einstein summation convention) is applied.

The first Piola-Kirchhoff stress tensor is deduced from the Helmholtz free energy function using

𝑷=∂ψ∂𝑭=∂ψ∂𝑪e:∂𝑪e∂𝑭.\boldsymbol{P}=\displaystyle{\frac{\partial\psi}{\partial\boldsymbol{F}}}=\displaystyle{\frac{\partial\psi}{\partial\boldsymbol{C}_{e}}}:\displaystyle{\frac{\partial\boldsymbol{C}_{e}}{\partial\boldsymbol{F}}}. (16)

Using Eq. 14, it can be deduced that

∂𝑪e∂𝑭=1ϑ2​∂𝑪∂𝑭.\displaystyle{\frac{\partial\boldsymbol{C}_{e}}{\partial\boldsymbol{F}}}=\displaystyle{\frac{1}{\vartheta^{2}}}\,\displaystyle{\frac{\partial\boldsymbol{C}}{\partial\boldsymbol{F}}}. (17)

Finally, the balance of linear momentum governs the quasi-static equilibrium of the arterial wall structure.

∇⋅𝑷+𝒃=𝟎\nabla\cdot\boldsymbol{P}+\boldsymbol{b}=\boldsymbol{0} (18)

3 Numerical methods

The Galerkin weak forms of the governing differential equations are arrived at.

∫Ω∂cP∂t​δ​cP​𝑑v\displaystyle\displaystyle{\int_{\Omega}\frac{\partial c_{{}_{P}}}{\partial t}}\,\delta c_{{}_{P}}\,dv =\displaystyle= ∫Γ(JP⋅𝒏)​δ​cP​𝑑a−∫Ω(∇δ​cP⋅DP​∇cP)​𝑑v\displaystyle\displaystyle{\int_{\Gamma}(J_{{}_{P}}\cdot\boldsymbol{n})\,\delta c_{{}_{P}}\,da}-\displaystyle{\int_{\Omega}(\nabla\delta c_{{}_{P}}\cdot D_{{}_{P}}\nabla c_{{}_{P}}})\,dv (19)
−\displaystyle- ∫Ωα​ρS​cP​δ​cP​𝑑v\displaystyle\displaystyle{\int_{\Omega}\alpha\,\rho_{{}_{S}}\,c_{{}_{P}}}\,\delta c_{{}_{P}}\,dv
∫Ω∂ρE∂t​δ​ρE​𝑑v\displaystyle\displaystyle{\int_{\Omega}\frac{\partial\rho_{{}_{E}}}{\partial t}}\,\delta\rho_{{}_{E}}\,dv =\displaystyle= ∫Ωβ​ρS​(1−ρEρE,t​h)​δ​ρE​𝑑v−∫Ωγ​cP​ρE​δ​ρE​𝑑v\displaystyle\displaystyle{\int_{\Omega}\,\beta\,\rho_{{}_{S}}\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)\,\delta\rho_{{}_{E}}}\,dv-\displaystyle{\int_{\Omega}\,\gamma c_{{}_{P}}\rho_{{}_{E}}\,\delta\rho_{{}_{E}}}\,dv (20)
∫Ω∂ρS∂t​δ​ρS​𝑑v\displaystyle\displaystyle{\int_{\Omega}\frac{\partial\rho_{{}_{S}}}{\partial t}}\,\delta\rho_{{}_{S}}\,dv =\displaystyle= ∫Γ(JS⋅𝒏)​δ​ρS​𝑑a−∫Ω(∇δ​ρS⋅χ​cP​(1−ρEρE,t​h)​ρS​∇ρE)​𝑑v\displaystyle\displaystyle{\int_{\Gamma}(J_{{}_{S}}\cdot\boldsymbol{n})\,\delta\rho_{{}_{S}}\,da}-\displaystyle{\int_{\Omega}\left(\nabla\,\delta\rho_{{}_{S}}\cdot\chi\,c_{{}_{P}}\,\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)\,\rho_{{}_{S}}\,\nabla\rho_{{}_{E}}\right)}\,dv (21)
+\displaystyle+ ∫Ωκ​cP​ρS​(1−ρEρE,t​h)​δ​ρS​𝑑v\displaystyle\displaystyle{\int_{\Omega}\kappa\,c_{{}_{P}}\,\rho_{{}_{S}}\left(1-\displaystyle{\frac{\rho_{{}_{E}}}{\rho_{{}_{E,th}}}}\right)}\,\delta\rho_{{}_{S}}\,dv
∫Ω0𝑷:δ​𝑭​d​V\displaystyle\displaystyle{\int_{\Omega_{0}}}\,\boldsymbol{P}:\delta\boldsymbol{F}\,dV =\displaystyle= ∫Γ0(𝑻⋅δ​𝒖)​𝑑A\displaystyle\displaystyle{\int_{\Gamma_{0}}}\,(\boldsymbol{T}\cdot\delta\boldsymbol{u})\,dA (22)

Here, JPJ_{{}_{P}} and JSJ_{{}_{S}} denote the PDGF and SMC fluxes across the boundary, 𝒏\boldsymbol{n} the normal on the boundary and 𝑻\boldsymbol{T} the traction on the boundary. Body forces are considered absent.

The transport equations are temporally discretised using the semi-implicit Backward-Euler method. Semi-implicitness is attributed to the fact that in each of the weak forms, variables other than the ones discretised in time are carried over from the previous time step. Spatial discretisation is performed using linear finite elements. A decoupled set of linear equations for the time-step updation of transport variables are hence obtained.

[𝑴+Δ​t​𝑳+Δ​t​𝑷]​𝒄Pn+1\displaystyle\left[\boldsymbol{M}+\Delta t\,\boldsymbol{L}+\Delta t\,\boldsymbol{P}\right]\boldsymbol{c}_{{}_{P}}^{n+1} =\displaystyle= 𝑴​𝒄Pn\displaystyle\boldsymbol{M}\boldsymbol{c}_{{}_{P}}^{n} (23)
[𝑴+Δ​t​𝑻]​𝝆En+1\displaystyle\left[\boldsymbol{M}+\Delta t\,\boldsymbol{T}\right]\boldsymbol{\rho}_{{}_{E}}^{n+1} =\displaystyle= 𝑴​𝝆En+𝑹\displaystyle\boldsymbol{M}\boldsymbol{\rho}_{{}_{E}}^{n}+\boldsymbol{R} (24)
[𝑴+Δ​t​𝑲−Δ​t​𝑸]​𝝆Sn+1\displaystyle\left[\boldsymbol{M}+\Delta t\,\boldsymbol{K}-\Delta t\,\boldsymbol{Q}\right]\boldsymbol{\rho}_{{}_{S}}^{n+1} =\displaystyle= 𝑴​𝝆Sn\displaystyle\boldsymbol{M}\boldsymbol{\rho}_{{}_{S}}^{n} (25)

Definitions of matrices in the equations above can be found in the appendix. The diffusion and advection operators within the finite element setting, when discretised with linear finite elements, display instabilities across steep gradients in the transport quantities. This leads to unphysical negative concentrations and densities at the diffusion and advection fronts. To treat this, artificial discrete diffusive fluxes are introduced, which result in over smoothened solution profiles. Hence antidiffusive fluxes, limited using the Zalesak limiter, are utilised to preserve local extrema kuzmin2002 ; strehl2013 .

The weak form of the balance of linear momentum in Eq. 22 is spatially discretised using the same linear finite elements as defined for the transport problems, and is solved iteratively using the Newton-Raphson method. The residual and its linearization, prior to discretisation, required for the iterative updates are determined as follows.

gu\displaystyle g_{{}_{u}} =\displaystyle= ∫Ω0𝑷:δ​𝑭​d​V−∫Γ0(𝑻⋅δ​𝒖)​𝑑A\displaystyle\displaystyle{\int_{\Omega_{0}}}\,\boldsymbol{P}:\delta\boldsymbol{F}\,dV-\displaystyle{\int_{\Gamma_{0}}}\,(\boldsymbol{T}\cdot\delta\boldsymbol{u})\,dA (26)
Δ​gu\displaystyle\Delta g_{{}_{u}} =\displaystyle= ∫Ω0δ𝑭:(𝒜:Δ𝑭)dV\displaystyle\displaystyle{\int_{\Omega_{0}}}\delta\boldsymbol{F}:(\mathcal{A}:\Delta\boldsymbol{F})\,dV (27)

The fourth order tensor 𝒜\mathcal{A} is computed using

𝒜=∂𝑷∂𝑭+∂𝑷∂ϑ⊗∂ϑ∂𝑭.\mathcal{A}=\displaystyle{\frac{\partial\boldsymbol{P}}{\partial\boldsymbol{F}}+\frac{\partial\boldsymbol{P}}{\partial\vartheta}\otimes\frac{\partial\vartheta}{\partial\boldsymbol{F}}}. (28)

Coupling

The transport problems are solved in the current configuration Ω\Omega while the balance of linear momentum is solved in the reference configuration Ω0\Omega_{0}. A staggered solution strategy is carried out to achieve coupling. The transport problem provides the SMC density to the structural problem while the structural problem updates the current configuration Ω\Omega based on the SMC density. An incremental form of Eq. 10 is applied in this regard as follows:

ϑ=ϑ−​(1+JJ−​(ρSρS,h−1))1/d,\vartheta=\vartheta^{-}\left(1+\displaystyle{\frac{J}{J^{-}}\,\left(\frac{\rho_{{}_{S}}}{\rho_{{}_{S,h}}}-1\right)}\right)^{1/d}, (29)

wherein the quantities ()−(\,\,)^{-} refer to those from the previous time step. At the end of solution of the structural problem, the transport quantities ϕ\phi are updated at the Gauss-points using ϕ∗=(J−/J)​ϕ\phi^{*}=(\,J^{-}/J)\,\phi and extrapolated to the nodes. The staggered solution strategy provides flexibility in case different spatial discretisations are absolutely necessary for the systems being coupled.

4 Numerical investigation

A longitudinal axisymmetric section of an artery is considered fixed on either ends. A Gaussian distribution of PDGF concentration consisting of three distinct peaks is prescribed in the domain to mimic a PCI injury in the arterial wall. ECM and SMC densities (ρE,0\rho_{{}_{E,0}} and ρS,0\rho_{{}_{S,0}}) are uniformly prescribed to be those of a healthy artery. Zero flux boundary conditions are prescribed for all the transport quantities on all boundaries. Some of the parameters are deduced from literature and the others tuned to predict the physics qualitatively.

Parameter Value [units]
\svhline DPD_{{}_{P}} 0.010.01 [mm2/day]
α\alpha 1.0×10−131.0\times 10^{-13} [mm3/cell/day]
β\beta 5.0×10−85.0\times 10^{-8} [mol/cell/day]
γ\gamma 5.0×10175.0\times 10^{17} [mm3/mol/day]
χ\chi 1.0×10191.0\times 10^{19} [mm5/cell/day]
κ\kappa 1.0×10−21.0\times 10^{-2} [mm3/mol/cell/day]
ρE,0\rho_{{}_{E,0}} 7.0×10−97.0\times 10^{-9} [mol/mm3]
ρE,t​h\rho_{{}_{E,th}} 1.1×ρE,01.1\times\rho_{{}_{E,0}}
ρS,0\rho_{{}_{S,0}} 3.16×1063.16\times 10^{6} [cells/mm3]
μa\mu^{a} 0.020.02 [M Pa]
λa\lambda^{a} 1010 [M Pa]
k1k_{1} 0.1120.112 [M Pa]
k2k_{2} 20.6120.61 [-]
κa\kappa^{a} 0.10.1 [-]
θia\theta^{a}_{i} 4141 [∘]
Table 2: Parameter set
Refer to caption
Figure 1: Problem setup

Transport simulation is run for a time period of 1 day to perform the mesh and time step convergence analyses since this dictates the spatial and temporal discretisations. The SMC density plots along section A-A are used for the analyses.

Refer to caption
Figure 2: mesh convergence (A-A)
Refer to caption
Figure 3: time-step convergence (A-A)

Using the converged discretisations, the coupled simulation is run for a period of 1 day and the results examined. It is expected that the growth ceases when the ECM is healed and PDGF is completely internalized. The evolutions of quantities at a Gauss-point of interest reflect this as shown in Fig. 4. Contraction is observed in some regions of the arterial wall which can be attributed to the migration of SMC away from these regions.

Refer to caption
Figure 4: deformed arterial wall and evolution of quantities of interest

5 Conclusion and outlook

The developed model is able to successfully capture the growth induced in the arterial wall due to mechanisms brought about by PDGF. It is to be noted that the model does not take into account the loss and healing processes of the endothelium. In reality, the restenotic process does not cease until the endothelium is completely healed and there is no further influx of PDGF. Wall shear stress dependent influx of PDGF can be prescribed as a boundary condition to achieve the aforementioned effect. The transforming growth factor (TGF)-β\beta is another significant constituent responsible for modulation of growth in vasculatures. It considerably affects the proliferative response of SMCs and hence needs to be considered in the framework developed in this work. Furthermore, it remains to be seen how sensitive the model is to the set of model parameters and which of those can be determined experimentally.

Future work shall therefore entail consideration of additional wall constituents that influence the restenotic process, physiological experimentation, parameter sensistivity evaluation, and finally arriving at a high-fidelity in-silico model which aides in examining the efficacy of drugs intended to be used on drug-eluting stents.

Acknowledgements.
The financial support of DFG for the projects ’Drug-eluting coronary stents in stenosed arteries: medical investigation and computational modelling’ (RE 1057 / 44-1) and ’Modelling of Structure and Fluid-Structure Interaction during Tissue Maturation in Biohybrid Heart Valves’ (RE 1057 / 45-1) is gratefully acknowledged.

Appendix

Element matrices in the linearised decoupled system of transport equations are listed here. 𝑵\boldsymbol{N} refers to the shape function vector associated with linear finite elements.

𝑴e=∫Ωe𝑵​𝑵T​𝑑v𝑳e=∫ΩeDP​∇𝑵​∇𝑵T​d​v𝑷e=∫Ωeα​ρSn​𝑵​𝑵T​𝑑v𝑻e=∫Ωe(β​ρSnρE,t​h+γ​cPn)​𝑵​𝑵T​𝑑v𝑲e=∫Ωe[χ​cPn​(1−ρEnρE,t​h)​∇𝑵​∇ρ1n​𝑵T]​𝑑v𝑸e=∫Ωeκ​cPn​(1−ρEnρE,t​h)​𝑵​𝑵T​𝑑v𝑹e=∫Ωeβ​ρSn​𝑵​𝑑v\displaystyle\begin{array}[]{r@{\mskip\thickmuskip}l}\boldsymbol{M}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\,\boldsymbol{N}\,\boldsymbol{N}^{T}\,dv\\ \boldsymbol{L}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\,D_{{}_{P}}\nabla\boldsymbol{N}\,\nabla\boldsymbol{N}^{T}\,dv\\ \boldsymbol{P}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\,\alpha\rho_{{}_{S}}^{n}\boldsymbol{N}\,\boldsymbol{N}^{T}\,dv\end{array}\quad\begin{array}[]{r@{\mskip\thickmuskip}l}\boldsymbol{T}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\left(\beta\,\displaystyle{\frac{\rho_{{}_{S}}^{n}}{\rho_{{}_{E,th}}}}+\gamma\,c_{{}_{P}}^{n}\right)\boldsymbol{N}\,\boldsymbol{N}^{T}\,dv\\ \boldsymbol{K}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\left[\chi\,c_{{}_{P}}^{n}\left(1-\displaystyle{\frac{\rho_{{}_{E}}^{n}}{\rho_{{}_{E,th}}}}\right)\nabla\boldsymbol{N}\,\nabla\rho_{{}_{1}}^{n}\boldsymbol{N}^{T}\right]\,dv\\ \boldsymbol{Q}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\,\kappa\,c_{{}_{P}}^{n}\left(1-\displaystyle{\frac{\rho_{{}_{E}}^{n}}{\rho_{{}_{E,th}}}}\right)\boldsymbol{N}\,\boldsymbol{N}^{T}\,dv\\ \boldsymbol{R}^{e}\mskip 5.0mu plus 5.0mu&=\displaystyle{\int_{\Omega_{e}}}\beta\,\rho_{{}_{S}}^{n}\,\boldsymbol{N}\,dv\end{array}

References

  • (1) Budu-Grajdeanu P, Schugart R C, Friedman A, Valentine C, Agarwal A K, Rovin B H et al (2008) A mathematical model of venous neointimal hyperplasia formation. Theor. Biol. Medical Model. 5:2
  • (2) Cyron C J, Aydin R C, Humphrey J D et al (2016) A homogenized constrained mixture (and mechanical analog) model for growth and remodeling of soft tissue. Biomech. Model Mechanobiol. 15:1389–1403
  • (3) Escuer J, Martínez, McGinty S, Peña et al (2019) Mathematical modelling of the restenosis process after stent implantation. J. R. Soc. Interface 16:20190313.
  • (4) Gasser T C, Ogden R W, Holzapfel G A et al (2006) Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface 3:15–35
  • (5) Himpel G, Kuhl E, Menzel A, Steinmann P et al (2005) Computational modelling of isotropic multiplicative growth. Comput Model Eng Sci 8(2):119–134
  • (6) Holzapfel G A, Gasser T C, Ogden R W et al (2000) A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models. Journal of Elasticity 61:1–48
  • (7) Humphrey J D, Rajagopal K R et al (2002) A constrained mixture model for growth and remodeling of soft tissues. Math. Models Methods Appl. Sci. 12(3): 407–430
  • (8) Keller E, Segel L A et al (1971) Model for Chemotaxis. J. theor. Biol. 30:225–234
  • (9) Koskinas K C, Chatzizisis S Y, Antoniadis A P, Giannoglou G D et al (2012) Role of Endothelial Shear Stress in Stent Restenosis and Thrombosis Pathophysiologic. Journal of the American College of Cardiology 59(15):1337–1349
  • (10) Koyama N, Hart C E, Clowes A W et al (1994) Different Functions of the Platelet-Derived Growth Factor-α\alpha and -β\beta Receptors for the Migration and Proliferation of Cultured Baboon Smooth Muscle Cells. Circulation Research 75:682–691
  • (11) Kuzmin D, Turek S et al (2002) Flux Correction Tools for Finite Elements. J. Comput. Phys. 175:525–558
  • (12) Li S, Lei L, Hu Y, Zhang Y, Zhao S, Zhang J et al (2019), A fully coupled framework for in silico investigation of in-stent restenosis. Comput. Methods Biomech. Biomed. Engin. 22(2):217–228
  • (13) Nolan D R, Gower A L, Destrade M, Ogden R W, McGarry J P et al (2014) J. Mech. Behav. Biomed. Mater. 39:48–60
  • (14) Rodriguez E K, Hoger A, McCulloch A D et al (1994) Stress-dependent finite growth in soft elastic tissues. J. Biomchanics 21(4):455–467
  • (15) Strehl R, Sokolov A, Kuzmin D, Horstmann D, Turek S et al (2013) A positivity-preserving finite element method for chemotaxis problems in 3D. J. Comput. Appl. Math. 239:290–303
  • (16) Thon M P, Hemmler A, Glinzer A, Mayr M, Wildgruber M, Zernecke-Madsen A, Gee M W et al (2018) A multiphysics approach for modeling early atherosclerosis. Biomech. Model Mechanobiol. 17:617–644
  • (17) Zahedmanesh H, van Oosterwyck H, Lally C et al (2012), A multi-scale mechanobiological model of in-stent restenosis: deciphering the role of matrix metalloproteinase and extracellular matrix changes. Comput. Methods Biomech. Biomed. Engin. 17(8):813–828