A Mechanobiological model for damage-induced growth in arterial tissue with application to in-stent restenosis

Accepted Manuscript

A Mechanobiological Model for Damage-induced Growth in Arterial Tissue with Application to In-stent Restenosis B. Fereidoonnezhad, R. Naghdabadi, S. Sohrabpour, G.A. Holzapfel PII: DOI: Reference:

S0022-5096(16)30824-9 10.1016/j.jmps.2017.01.016 MPS 3052

To appear in:

Journal of the Mechanics and Physics of Solids

Received date: Revised date: Accepted date:

12 November 2016 12 January 2017 18 January 2017

Please cite this article as: B. Fereidoonnezhad, R. Naghdabadi, S. Sohrabpour, G.A. Holzapfel, A Mechanobiological Model for Damage-induced Growth in Arterial Tissue with Application to In-stent Restenosis, Journal of the Mechanics and Physics of Solids (2017), doi: 10.1016/j.jmps.2017.01.016

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A Mechanobiological Model for Damage-induced Growth in Arterial Tissue with Application to In-stent Restenosis B. Fereidoonnezhada , R. Naghdabadia,b , S. Sohrabpoura , G.A. Holzapfelc,d, 1 a

c

Institute for Nano-Science and Technology, Sharif University of Technology, Tehran, Iran

CR IP T

b

Department of Mechanical Engineering, Sharif University of Technology, Tehran, Iran

Institute of Biomechanics, Graz University of Technology, Stremayrgasse 16-II, 8010 Graz, Austria

d Faculty

of Engineering Science and Technology, Norwegian University of Science and Technology 7491 Trondheim, Norway

To Appear in J. Mech. Phys. Solids

Abstract.

AN US

January 19, 2017

In-stent restenosis (ISR) is one of the main drawbacks of stent implementation

which limits the long-term success of the procedure. Morphological changes occurring within the arterial wall due to stent-induced mechanical injury are a major cause for activation of

M

vascular smooth muscle cells (VSMCs), and the subsequent development of ISR. Considering the theory of volumetric mass growth and adopting a multiplicative decomposition of the

ED

deformation gradient into an elastic part and a growth part, we present a mechanobiological model for ISR. An evolution equation is developed for mass growth of the neointima, in which

PT

the activation of VSMCs due to stent-induced damage (injury) and the proliferation rate of the activated cells are considered. By introducing the mass evolution into the mass balance

CE

equation, we obtain the evolution of the growth tensor over time. The model is implemented in a finite element code and the procedure of angioplasty is simulated, whereby the features

AC

of the proposed growth model are illustrated. Keywords:

Mechanobiological model, finite element simulation, stent deployment, damage-

induced growth, in-stent restenosis 1

To whom correspondence should be addressed. Email address: [email protected]

1

ACCEPTED MANUSCRIPT

1. Introduction Despite the great developments in medical and surgical treatment of cardiovascular disease (CVD), it has remained the major cause of death worldwide [1]. Atherosclerosis is one (major) type of CVD in which fibrous and fatty materials, called plaque, build up inside the artery

CR IP T

and cause partial or total occlusion of the artery (see, e.g., [2]). One way to restore the physiological blood flow of an occluded artery involves the deformation of plaque using intravascular balloon angioplasty with or without stenting. This clinical intervention is accompanied with supra-physiological loading of the arterial wall causing tissue damage (injury) (see [3] for more details). It may also lead to tissue in-growth and re-blockage, termed in-stent restenosis

AN US

(ISR) [4], which is one of the major complications associated with stent implementation, often limiting the long-term success of the procedure; usually ISR occurs within six months after stent deployment (see, e.g., [5]). The mechanism of restenosis after angioplasty is a combination of elastic recoil, arterial vessel remodeling, and neointimal hyperplasia. Late lumen loss in stented segments of the artery is often the result of neointimal hyperplasia [6].

M

Despite significant improvements on the stent design, up to 10% of the stenting procedures require a repetition due to ISR, a rate that increases up to 50% in some high-risk lesions [4].

ED

Drug-eluting stents (DESs), which are coated with anti-mitogenic or anti-proliferative drugs to prevent smooth muscle cell proliferation and neointimal growth, decrease the rate of ISR in

PT

comparison with bare metal stents in coronary arteries [7]. However, late in-stent thrombosis and lack of re-endothelialization are the major drawbacks associated with the use of DESs [8].

CE

A successful design of stents requires a mechanobiological model that can predict ISR induced as a function of mechanical changes in the arterial tissue after stent deployment. In order to

AC

develop such a predictive model it is important to consider the stent–induced changes on the cellular level, that is excessive migration and proliferation of vascular smooth muscle cells (VSMCs) [9]. Given that the degree of vessel injury induced by stent insersion is a stimulus for ISR [10–13] it is a key design aspect for stents. In this context, computational modeling techniques are efficient and powerful to optimize stent design parameters in order to minimize stent–induced arterial stresses and injury. Several computational models of stent 2

ACCEPTED MANUSCRIPT

deployment procedures have been developed in recent years with the goal of reducing stentinduced arterial stresses [14–20]. The influence of several mechanical parameters such as the stent strut thickness [15, 21, 22] and the plaque composition [23, 24] have been studied. Balloon expandable stents have also been compared with self-expanding stents in terms of the level of stresses they induce within the arterial wall, and hence the risk of arterial injury using

CR IP T

finite element (FE) models [25].

Moreover, computational modeling is a powerful approach to capture the biological response of an implant using mechanobiological models which relate the mechanical environmental changes to the growth and remodeling of vascular cells [4, 26]. The approaches utilizing discrete methods such as cellular automata [27, 28] and agent-based models [29, 30] have

AN US

been employed for modeling ISR [8, 31, 32]. The approach adopted in these studies is mainly based on discrete methods, whereas a continuum method using the FE method may also be considered as valuable. FE methods have shown to be robust for the quantification of arterial stresses, and they have been successfully utilized for patient-specific modeling of the stent–

M

artery interaction [20, 33].

We need to know more about ISR. A constitutive model based on the cellular mechanisms ac-

ED

tivated after stent deployment, can improve the prediction of ISR. In addition, the numerical implementation of such a mechanobiological model provides a powerful approach to predict late lumen loss due to ISR, and it also optimizes the stenting procedure in terms of devices

PT

(such as stent, catheter, etc.) and procedure parameters (such as inner balloon pressure, bal-

CE

loon inflation time, etc.).

The main goal of the present study is to propose a continuum and computational model for

AC

ISR that takes into account the cellular changes in the arterial wall caused by stent-induced injury. The underlying idea is the use of a multiplicative decomposition of the deformation gradient into an elastic and a growth part within the framework of volumetric finite growth. In this context, the characterization of the growth tensor evolution is one of the major challenges, where, by most of the available models, a phenomenological evolution for the growth tensor [34, 35]. Although the mathematical understanding of these evolution equations seems to be

3

ACCEPTED MANUSCRIPT

fairly intuitive, their clear biomechanical interpretation is still widely debated and unclear. Thus, it seems natural to seek for a deeper insight into the micromechanical phenomena that are related to tissue growth on the cellular level. Motivated by cellular changes after stent implementation (migration and proliferation of VSMCs), we propose an evolution equation for the mass of the neointima. Considering the effect of the damage intensity on the number of

CR IP T

activated VSMCs, and the proliferation rate of the activated VSMCs, the evolution equation has a clear biomechanical interpretation. The evolution for the growth tensor is then obtained by introducing the mass evolution into the balance equation of mass.

The paper is organized as follows: Section 2 explains the mathematical framework and introduces the continuum model for finite growth. The kinematics and mass balance equation are

AN US

presented, and thermodynamic consistency are investigated. Furthermore, the general constitutive model is specialized for ISR of arterial tissue in terms of an energy function, a density evolution and a growth tensor. In Section 3, the FE implementation and the related verification are presented. The stress and elasticity tensors are determined, and the solution algorithm is

M

explained. A parameter study of the model is conducted in Section 4. The procedure of angioplasty with subsequent restenosis is simulated, and the results are presented in Section 5.

ED

In Section 6 we draw some conclusions.

PT

2. Modeling of finite growth

In this section we illustrate the governing equations for finite growth which consist of: (i) the

CE

kinematic equations, which utilize the concept of an incompatible growth configuration, and the multiplicative decomposition of the deformation gradient; (ii) the mass balance equation,

AC

which is an important part of the growth process; (iii) the derivation of the reduced inequality showing thermodynamic consistency; (iv) the constitutive equations, which define the baseline elastic response in the intermediate configuration, supplemented by (v) an appropriate form for the growth tensor and its evolution equation to characterize the ISR phenomenon.

4

ACCEPTED MANUSCRIPT

2.1. Kinematics Within the framework of finite growth, the key kinematic assumption is based on the consideration of an intermediate configuration Bg , i.e. between the initial configuration B0 and the current configuration Bt , which naturally leads to the multiplicative decomposition of the

F = Fe Fg .

CR IP T

deformation gradient F into an elastic part Fe and a growth part Fg , [36]. Thus, (1)

We also note that the arterial layers in the initial configuration are characterized by two collagen fiber families, symmetrically distributed by an angle β with respect to the circumferential direction of the artery. Thus, the arterial layers respond orthotropically, and the families of

AN US

collagen fibers are characterized by the two direction vectors a10 and a20 with |a10 | = 1 and |a20 | = 1, in both the initial and the intermediate configuration (assuming isotropic growth).

The Jacobian of the deformation gradient F and its elastic and growth counterparts are denoted by J = det F, Je = det Fe , and Jg = det Fg , respectively, such that

M

J = Je Jg .

(2)

Next we introduce the right Cauchy–Green tensor C= FT F and the elastic right Cauchy–

ED

Green tensor, defined by Ce = FT e Fe , as

PT

−1 Ce = F−T g CFg .

(3)

CE

˙ −1 then takes on the following form The velocity gradient l= FF −1 l = F˙ e Fe −1 + Fe (F˙ g F−1 g )Fe .

(4)

The pull-back operation of the spatial velocity gradient (4) to the intermediate configuration

AC

Bg leads

L = F−1 e lFe = Le + Lg ,

˙ Le = F−1 e Fe ,

with

Lg = F˙ g F−1 g .

(5)

Moreover, from (4) and the property of the trace of a matrix product, the following additive decomposition is obtained ˙ −1 trl = tr(F˙ e F−1 e ) + tr(Fg Fg ). 5

(6)

ACCEPTED MANUSCRIPT

Subsequently, we also consider the multiplicative split Ce = Je 2/3 Ce ,

(7)

where Ce and Je represent the isochoric and volumetric part of the elastic deformation, respectively.

CR IP T

2.2. Mass balance equation

Unlike engineering materials, living biological tissues are known to adapt (grow/atrophy) in response to an environmental change, whereby the mass of the system is changed, and consequently the consideration of the mass balance equation is an important point for such

AN US

systems. If the infinitesimal mass dm0 of an infinitesimal volume element dV0 in the initial configuration is dm0 = ρ0 dV 0 , where ρ0 denotes the related density, the infinitesimal mass dm of the corresponding element in the intermediate and current configurations is dm = ρg dVg = ρdV ,

(8)

M

where ρg and dVg denote the density and the infinitesimal volume element in the intermediate configuration, and ρ and dV are the corresponding quantities in the current configuration. It

ED

is worth noting that the growth, and consequently the mass change of the biological tissue, occurs between the initial and the intermediate configuration, whereby the mass of an element in the current configuration is the same as the mass of the corresponding element in the

PT

intermediate configuration.

CE

Subsequently, it is convenient to designate the quantity ρJ as ρ0g , where J = dV /dV0 . Thus, ρ0g =

dm = ρJ. dV 0

(9)

AC

In addition, with eq. (8) it follows that ρ0g = ρg Jg ,

ρg = ρJe ,

(10)

where Jg = dVg /dV0 , and eq. (2) was used. Next, if rg and rg0 are time rates of the mass growth per unit current volume V and initial volume V0 , respectively, then d (dm) = rg dV = rg0 dV0 , dt 6

(11)

ACCEPTED MANUSCRIPT

where d/dt(•) stands for the material time derivative of (•). In addition, in view of (9)1 and (11)2 we can conclude that dρ0g = rg0 . dt

(12)

0

CR IP T

By integrating eq. (11)2 , the following expression is obtained for the mass balance, i.e. Z t rg0 dτ dV0 . (13) dm = dm0 + In view of (9)2 and (10)1 , the mass balance equation (13) may then be represented as Z t ρg Jg = ρ0 + rg0 dτ . 0

From (14) and (10) we further find that

AN US

d d (ρJ) = (ρg Jg ) = rg0 . dt dt

(14)

(15)

This yields the continuity equations for the densities ρ and ρg according to dρg + ρg tr(F˙ g F−1 g ) = rg Je . dt

dρ ˙ −1 ) = rg , + ρtr(FF dt

(16)

dρ dρg + ρg tr(F˙ e F−1 . e ) = dt dt

(17)

ED

Je

M

In view of (6), from eq. (16), it readily follows that

For the growth of soft biological tissues, we assume that the newly born tissue has the same density as the existing tissue (ρg = ρ0 = constant), an assumption which is frequently used

PT

in the literature, see, e.g., [37]. With this assumption we have ρJ = ρ0 Jg , ρ0 = ρJe , and the

CE

continuity equations (16) may be represented as dρ + ρtr(F˙ e F−1 e ) = 0, dt

ρtr(F˙ g F−1 g ) = rg .

(18)

dρ = 0, dt

(19)

AC

Moreover, due to incompressibility we have Je = 1,

tr(F˙ e F−1 e ) = 0,

and

ρ = ρ0 = constant.

Consequently, the continuity equation for an incompressible tissue, which grows in a density preserving manner, may be represented as tr(F˙ g F−1 g ) = 7

rg0 rg = 0. ρg ρg

(20)

ACCEPTED MANUSCRIPT

2.3. Thermodynamic consistency For a system with varying mass we now consider the Clausius–Duhem inequality for an isothermal process. Hence, 1 ˙ − ρ0 T S0 ≥ 0, S : C˙ − ρ0g Ψ g 2

(21)

CR IP T

where S is the second Piola–Kirchhoff stress tensor, T is the absolute temperature and the ˙ The extra entropy source S0 in the first derivative of Ψ with respect to time is denoted by Ψ. dissipation inequality, which takes into account the increase of entropy due to non-mechanical components, is necessary to satisfy thermodynamic consistency in the context of growth [34]. Inequality (21) must be fulfilled for arbitrary thermodynamic processes to guarantee the non-

AN US

negativeness of the internal dissipation (required by the second law of thermodynamics).

Inequality (21) is general and does not include any indication of the variables on which Ψ depends. At this point, in order to obtain an explicit constitutive law, we need to specify such variables. We assume an elastic response at any stage of the growth and that the stress depends

M

on the growth. Damage-induced inelastic phenomena such as stress softening and permanent deformation, as proposed in [38], can be further included in a rather straightforward manner.

ED

The free energy per unit current mass of the biological tissue can be expressed in terms of the elastic right Cauchy–Green tensor Ce , the structural tensors A10 = a10 ⊗ a10 and A20 = a20 ⊗ a20

Ψ = Ψ(Ce , A10 , A20 , ρ0g ).

(22)

CE

PT

and the density ρ0g as

By differentiating the free-energy function (22) with respect to time, we have

AC

∂ρ0 ∂ρ0 ˙ = ∂Ψ : ∂Ce : C˙ + ∂Ψ : ∂Ce : F˙ g + ∂Ψ g : F˙ g + ∂Ψ g ρ˙ g Ψ ∂Ce ∂C ∂Ce ∂Fg ∂ρ0g ∂Fg ∂ρ0g ∂ρg ∂Ψ ∂Ψ ˙ 2 1 + 1 : A˙ 0 + : A0 , ∂A0 ∂A20

(23)

in which C, Fg , ρg , A10 and A20 are considered to be the independent variables, and Ce and ρ0g are related to the independent variables by eqs. (3) and (10)1 , respectively. Then, in view of 8

ACCEPTED MANUSCRIPT

eqs. (3) and (10)1 , eq. (23) can be represented as (for a detailed derivation see Appendix A) ˙ = F−1 ∂Ψ F−T : C˙ − 2Ce ∂Ψ : F˙ g F−1 + ρ0 ∂Ψ I : F˙ g F−1 Ψ g g g g ∂Ce g ∂Ce ∂ρ0g ∂Ψ ∂Ψ ˙ 2 ∂Ψ ˙ 1 +Jg 0 ρ˙ g + : A0 . 1 : A0 + ∂ρg ∂A0 ∂A20

(24)

Since we are considering no fiber reorientation, the structural tensors A10 and A20 change only

CR IP T

1 2 through F, and no separate evolution law is required. Thus A˙ 0 = A˙ 0 = 0. From eq. (24) and

AN US

the inequality (21) we then find that     ∂Ψ 1 0 −1 ∂Ψ −T 0 0 ∂Ψ ˙ S − ρ g Fg F : C + ρg 2Ce − ρg 0 I : F˙ g F−1 g 2 ∂Ce g ∂Ce ∂ρg ∂Ψ −ρ0g Jg 0 ρ˙ g − ρ0g T S0 ≥ 0. ∂ρg

(25)

Thus, with the standard arguments of continuum thermodynamics we obtain the following expression for the second Piola–Kirchhoff stress tensor, i.e. Se = 2ρ0g

∂Ψ . ∂Ce

(26)

M

−T S = F−1 g Se Fg ,

Moreover, we note that for the density preserving growth, which seems reasonable in the

ED

context of arterial growth, the density in the intermediate configuration is constant (ρ˙ g = 0), and with the constitutive relation (26)2 and with (5)3 the reduced dissipation inequality (25) gives

PT

  0 2 ∂Ψ Me − (ρg ) I : Lg − ρ0g T S0 ≥ 0, ∂ρ0g

(27)

CE

where the Mandel stress tensor Me = Ce Se of the intermediate configuration has been introduced. We will return to this inequality in Section 2.6, after the specification of the free-energy

AC

function Ψ, the growth tensor Fg and its evolution. 2.4. Specific form of the free-energy function Various forms of the free-energy function for arterial tissues were proposed in the literature. The model in [39] is one which is frequently used and implemented in several commercial

9

ACCEPTED MANUSCRIPT

software codes such as

ABAQUS

[40] and

FEAP

[41]. Employing this particular model, we

ˆ i.e. consider the following structure of the free energy per unit initial volume, say Ψ, ˆ = Ψ(C ˆ e , A1 , A2 ) = Ψ ˆ vol (Je ) + Ψ ˆ iso (Ce , A1 , A2 ), ρ0g Ψ = Ψ 0 0 0 0

(28)

ˆ vol (Je ) and Ψ ˆ iso (Ce , A1 , A2 ) describe the volumetric and the isochoric elastic response where Ψ 0 0

CR IP T

of the biological tissue. It is noted that the additive decomposition of the free-energy function into isochoric and volumetric parts has some limitations when applied to materials which are not incompressible (see, e.g., [42, 43]), which is not the case for arterial tissues.

ˆ iso = Ψ ˆ iso + Ψ m

X

ˆ iso , Ψ f,i

i=4,6


ˆ iso ¯ Ψ m = µ I1e − 3 ,

AN US

The isochoric part of the free-energy function (28) may be represented as


 k1  ˆ iso Ψ exp k2 E¯i2 − 1 , f,i = 2k2

(29)

i = 4, 6,

(30)

where i = 4, 6 relates to the two fiber families, µ, k1 , and k2 are material parameters, and I¯1e = trCe is the first invariant of Ce . In (30)2 E¯i = κI¯1e + (1 − 3κ)I¯ie − 1, i = 4, 6, denotes

M

the Green–Lagrange strain–like quantity, where I¯4e = a10 · Ce a10 and I¯6e = a20 · Ce a20 are the

modified pseudo-invariants, and κ ∈ [0, 1/3] characterizes the dispersion parameter. Note that

ED

the explicit expression in (30)2 can be replaced by other functions as, e.g., provided in [44]

which considers a non-symmetric dispersion of the collagen fibers better reflecting structural

(31)

CE

PT

ˆ vol according to data [45]. In addition, we consider the convex form of Ψ  2  1 J − 1 e vol ˆ (Je ) = Ψ − ln Je , D1 2

AC

which we have numerically used as a penalty function, where 1/D1 is a penalty parameter. 2.5. Micromechanically motivated evolution for the mass In this section we propose a mathematical model for mass growth of the intima triggered on the cellular level after angioplasty with stent deployment. The process of angioplasty with stent deployment denudes the endothelial layer and may stretch, injure or even rupture parts of tissue layers. This arterial injury activates platelet aggregation, triggers the migration of 10

ACCEPTED MANUSCRIPT

the VSMCs from the medial wall towards the intima where they start to proliferate and form a newly born material, called neointima. Under physiological conditions, VSMCs tend to remain in a contractile quiescent state (G0 of the cell cycle). However, after arterial injury, some of the VSMCs change their phenotypes and enter the G1 growth phase so that the cells become less adhesive to its neighborhood in the later phase of the G1 cell cycle, provoking

CR IP T

cell migration [9, 46].

Based on these cellular changes, we assume that the time-dependent rate of the mass growth of the neointima is dependent on the extent of injury to the vascular wall induced by the stent, and the current number of proliferating cells in the population. As an evolution equation

expression rg0 =

AN US

for the current mass per unit initial volume ρ0g of the neointima we propose the following dρ0g = f g (ρ0g , D, t), dt

(32)

where D is a measure of damage, t stands for time and f g denotes a growth function. Next,

M

we decompose f g as

f g (ρ0g , D, t) = f1g (ρ0g , t)f2g (D),

(33)

ED

where function f1g represents the time variation of mass due to the change in the cell proliferation rate, and f2g takes into account the effect of damage on the number of activated VSMCs.

PT

In other words, the number of VSMCs which migrate to the arterial inner surface are directly related to damage/injury induced in the tissue, and the function f2g takes into account such a

CE

dependency.

Moreover, the function f1g (ρ0g , t) shows how the migrated VSMCs proliferate or die. The

AC

mathematical model of proliferative cell kinetics can be adopted to represent this functionality. Here, in analogy with the cell proliferating model introduced in [47], for f1g we propose the following expression f1g (ρ0g , t) = kρ0g t exp(−αt),

(34)

where α and k are material parameters, and ρ0g is the time–dependent value of mass per unit initial volume of the cells in the population according to (9). This expression contains two 11

ACCEPTED MANUSCRIPT

parts: (i) the first term kρ0g t, which provides cellular growth within the population proportional to the mass of the cells already present; (ii) the second (exponential) term accounts for the decrease in proliferation within the population over time both due to cell death and daughter cells that cease dividing [47]. We highlight that the first term in (34) is an ascending function, while the second term is a descending function which implies two cellular mechanisms (cell

CR IP T

proliferation and death). For small values of t the first term is dominant and the growth rate increases up to a specific time, afterwards the second exponential term becomes dominant and the growth rate decreases and approaches to zero. This is consistent with clinical observations. Finally, the explicit form of the function f2g must be specified. First, a quantification of damage within the arterial tissue is required. One way to calculate damage is to use experimental data,

AN US

to establish a relationship between a measurable quantity (e.g., strain) or calculable quantity (e.g., stress) and damage in the arterial wall. In the absence of such experimental data, an ˆ iso alternative way is to choose the maximum isochoric free energy, say Ψ max , attained in the history of deformation as a (phenomenological) measure of damage, as discussed in [48] and

M

[38]. Following [38] we may define the damage D as   1 1 ˆ iso D = erf Ψ , r1 m1 max

(35)

ED

where erf(•) is the error function of (•), and r1 , m1 are damage parameters. Then for f2g we

PT

propose the following expression

f2g (D) = hD − Dth i,

(36)

where Dth is a threshold value for the growth to be started, and h i stands for the Macaulay

CE

brackets. It indicates that the growth occurs only when the damage in the arterial wall exceeds a certain threshold value.

AC

Finally, in view of (34) and (36), the growth function (33) takes on the form f g (ρ0g , D, t) = kρ0g t exp(−αt)hD − Dth i.

(37)

2.6. Specific form for the growth tensor and its evolution To complete the constitutive formulation, a specific form for the growth tensor Fg is presented. The growth of the anisotropic arterial tissue is expected to be anisotropic. However, to the 12

ACCEPTED MANUSCRIPT

authors’ knowledge to date no experimental data are available to characterize the anisotropic neointimal growth so that for simplicity we start with an isotropic growth tensor of the form Fg = λg I,

(38)

the velocity gradient (5)3 in the intermediate configuration is λ˙ g Lg = F˙ g Fg −1 = I. λg

CR IP T

where λg is the growth stretch and I is the second-order identity tensor. It readily follows that

(39)

In view of (32)2 and (37), the continuity equation (20)3 for an incompressible tissue can then be represented as

AN US

3λ˙ g = kt exp (−αt) hD − Dth i, λg

(40)

which gives the evolution of the growth stretch as

(41)

M

kt λ˙ g = exp(−αt)λg hD − Dth i. 3

Here we note that the damage parameter D in (41) is a measure of injury induced in the tissue

ED

during angioplasty and stenting, and it acts as a trigger for VSMC migration. It is calculated at the end of the procedure and is not a time-dependent quantity. By integrating eq. (41) by parts and by using the boundary condition λg = 1 at time t = 0, we obtain after a (lengthy

CE

PT

but) straightforward manipulation a closed form solution for λg , i.e.   k th [1 − (1 + αt) exp(−αt)]hD − D i , λg = exp 3α2

(42)

which provides the growth stretch λg over time for a growing restenotic lesion. This is illus-

AC

trated in Fig. 1(a) for arbitrary values of parameters (summarized in Table 1; we also used ˆ iso = 0.02 MPa), while the growth stretch rate λ˙ g is shown in Fig. 1(b). Ψ max

As can be seen from Fig. 1(b) the growth stretch rate λ˙ g increases up to a maximum value, and after that it decays until the growth stretch λg approaches its maximum value, i.e.   k max th λg = exp hD − D i . 3α2 13

(43)

ACCEPTED MANUSCRIPT

3.5

Growth stretch rate λ˙ g (1 · 103 )

1.05

λmax g

Growth stretch λg

1.04

1.03

1.02

1.01

5

10

15

20

25

30

35

40

45

2.5 2 1.5 1 0.5 0 0

50

5

10

15

20

25

30

35

40

45

50

CR IP T

1 0

3

Time (days)

Time (days)

(a)

(b)

Figure 1: Plots of (a) the growth stretch λg and (b) the growth stretch rate λ˙ g as a function of time (in days) for ˆ iso = 0.02 MPa. The growth stretch λg the parameters presented in Table 1, and a maximum free energy of Ψ max according to (43), while the growth stretch rate λ˙ g increases to a increases monotonically and approaches λmax g

AN US

maximum value and then decreases to zero.

Table 1: Material parameters for the proposed model used for verification and for a parameter study.

Parameter

Value

Unit

Hyperelastic parameters

k1 k2

MPa

0.56

MPa

16.21

−

51

◦

κ

0.18

−

ED

β

PT CE AC

0.0085

M

µ

Damage parameters

−

r1

1.59

m1

0.008

MPa

Dth

0.1

−

Growth parameters α

0.2

day−1

k

0.01

day−2

This non-monotonic rate of growth is consistent with clinical observation, see, e.g., [47]. Next let us find the time at which λ˙ g is a maximum. The time derivative of (41) results to   2  dλ˙ g k kt th th = exp(−αt)hD − D iλg exp(−αt)hD − D i + 1 − αt . dt 3 3 14

(44)

ACCEPTED MANUSCRIPT

Hence, we need to find the roots of (44). The first bracket is zero when D ≤ Dth or t

approaches infinity, where for both cases λ˙ g = 0. Thus, the time at which the growth velocity is a maximum can be obtained by setting the second bracket in (44) equal to zero, i.e. 1 − αt +


kt2 exp(−αt) D − Dth = 0, 3

(45)

discussed above.

CR IP T

where the Macaulay brackets are eliminated because (45) has to be solved for D > Dth , as

Solving the reduced dissipation inequality (27) for the extra entropy term we use Lg = λ˙ g I/λg , 0 2 ˆ i.e. eq. (39)2 , and (28)1 which leads to ∂Ψ/∂ρ0g = −Ψ/(ρ g ) . Thus we get

1 λ˙ g ˆ (trMe + 3Ψ). ρ0g T λg

AN US

S0 ≤

(46)

Remark: By using the multiplicative decomposition of the deformation gradient F we note that the intermediate configuration is not unique, in general. However, with the symmetry of Fg in (38), the rotation is lumped into the elastic part of the deformation gradient, and the

M

intermediate configuration is determined uniquely.

ED

3. Finite element implementation

PT

Due to the nonlinear governing equations required for finite growth, it is virtually impossible to find analytical solution for a boundary-value problem such as angioplasty with stenting.

CE

To solve the governing equations numerically, the FE method is employed to discretize the equations for finite growth in space. In this section, we first derive the expressions for the stress and elasticity tensors, then we describe the solution algorithm, and finally we provide a

AC

verification of the implemented model. 3.1. Stress tensor Once we have determined the growth stretch λg for a given deformation gradient F, we can successively determine the growth deformation gradient Fg from (38), the elastic deformation

15

ACCEPTED MANUSCRIPT

gradient Fe , the elastic right Cauchy–Green tensor Ce from (3) and the modified right Cauchy– Green tensor Ce from (7). Furthermore, it follows from (26)2 and (28) that Se = Svol e + Se ,

Svol e = 2

ˆ vol ∂Je ∂Ψ , ∂Je ∂Ce

Se = 2

ˆ iso ∂Ce ∂Ψ : , ∂Ce ∂Ce

(47)

where Svol e and Se are the volumetric and isochoric parts of the elastic second Piola–Kirchhoff

CR IP T

stress tensor Se , respectively. With the specific volumetric free energy (31), the volumetric part is further derived as −1 Svol e = pJe Ce ,

p=

1 (Je − Je−1 ), D1

(48)

where the relation ∂Je /∂Ce = Je C−1 e /2 has been used, and p denotes the hydrostatic pressure. is represented as Se =

m Se

+

X

f,i Se ,

AN US

In addition, with (29), the isochoric part of the elastic second Piola–Kirchhoff stress tensor Se

m Se

i=4,6

ˆ iso ∂Ψ =2 m , ∂Ce

f,i Se

ˆ iso ∂Ψ f,i . =2 ∂Ce

(49)

M

On use of the free-energy function (30), the first and second parts of the second Piola– Kirchhoff stress tensor in (49) are given by

ˆ iso ˜ m = 2 ∂ Ψm = 2µI, S e ∂Ce

ED

m ˜m, Se = Je −2/3 P : S e

PT

and

CE

f,i ˜ f,i , Se = Je −2/3 P : S e

AC

with

˜ f,i = 2 S e

ˆ iso ∂Ψ f,i = k1 E¯i exp(k2 E¯i2 ), ¯ ∂ Ei

ˆ iso ∂Ψ f,i ∂Ce

=2

ˆ iso ∂ E¯ ∂Ψ i f,i , ∂ E¯i ∂Ce

∂ E¯i = κI + (1 − 3κ)Ai0 , ∂Ce

(50)

(51)

(52)

where the fourth–order projection tensor P = I − C−1 e ⊗ Ce /3 has been introduced, while I

is the fourth–order unit tensor [49]. Once the elastic second Piola–Kirchhoff stress tensor Se

has been determined from equations (47)–(52), the total stress tensor S can then be obtained from (26)1 . Moreover, the Kirchhoff stress tensor τ is obtained by a push–forward operation of the second Piola–Kirchhoff stress tensor according to τ = FSFT . 16

ACCEPTED MANUSCRIPT

3.2. Elasticity tensor Consistent linearization of the second Piola–Kirchhoff stress tensor S with respect to the right Cauchy–Green tensor C gives the Lagrangian elasticity tensor C, which is essential for the application of the proposed (nonlinear) constitutive model within the FE method. Thus, with

Appendix B) C=2

CR IP T

(26)1 , (3) and the use of the index notation we obtain (for a more detailed derivation see

−T ∂(F−1 ∂S g Se Fg ) −1 −T −T =2 = (F−1 g ⊗Fg ) : Ce : (Fg ⊗Fg ), ∂C ∂C

with ∂Se iso Ce = 2 = Cvol e + Ce , ∂Ce

∂Svol =2 e , ∂Ce

Ciso e = 2

AN US

Cvol e

∂Se , ∂Ce

(53)

(54)

where the symbol ⊗ denotes a non–standard tensor product between two second-order tensors

according to (A⊗B)ijkl = Aik Bjl (see, e.g., [50]). calculated as [49]

The volumetric contribution is then

−1 −1 −1 Cvol ¯C−1 e = p e ⊗ Ce − 2pCe Ce ,

M

p¯ = p + Je

dp , dJe

(55)

ED

where the symbol has been introduced to denote the tensor product according to the rule −1 −1 −1 −1 −1 1 iso (C−1 e Ce )IJKL = 2 (Ce IK Ce JL + Ce IL Ce JK ). Moreover, with (49), the isochoric part Ce

PT

is obtained in the form of

m

Ciso e = Ce +

X

Ce ,

f,i Ce

∂S =2 e ∂Ce

f,i

(56)

i=4,6

AC

CE

where the two abbreviations

m

m Ce

f,i

∂S =2 e , ∂Ce

(57)

have been introduced. It is now required to provide explicit expression for the isochoric part Ciso e to be used in the FE implementation. Appendix C provides the related details of the derivations. The related Eulerian elasticity tensor c is then obtained by a standard push–forward operation

of C according to (C)ijkl = FiI FjJ FkK FlL (C)IJKL . For the isotropic growth tensor (38) and 17

ACCEPTED MANUSCRIPT

the multiplicative decomposition (1), (C)ijkl can be represented in the following modified form (C)ijkl = λ4g Fe iI Fe jJ Fe kK Fe lL (C)IJKL .

(58)

On using the isotropic growth tensor (38) and eq. (53)3 , we may deduce from (58) after some straightforward recasting in index notation the components of the Eulerian elasticity tensor C

CR IP T

in its final form as

(59)

(C)ijkl = Fe iI Fe jJ Fe kK Fe lL (Ce )IJKL .

In other words, for the proposed growth model, the Eulerian elasticity tensor

C

is obtained

by a push-forward of Ce via the elastic deformation gradient Fe . Note, however, that an

AN US

implementation into the commercial nonlinear FE software ABAQUS / STANDARD [40] requires a constitutive relationship in terms of the co-rotational elasticity tensor Λ, [51]. The link between the components of Λ and the components of C is given by

(60)

M

J(Λ)ijkl = J(C)ijkl + (τil δjk + τjk δil + τik δjl + τjl δik )/2,

3.3. Solution algorithm

ED

where δij is the Kronecker delta; for more details see, e.g., [38].

Table 2 illustrates the algorithmic flowchart for damage-induced isotropic growth. Remark-

PT

ably, throughout the entire algorithm, we never have to compute or store the coefficients of the growth tensor Fg . Instead, the algorithm manifests itself in a simple scalar scaling with

CE

λ−1 g , and is therefore straightforward and easy to implement. Moreover, we highlight that a closed form solution is presented for the growth stretch λg in eq. (42), which eliminates the

AC

local Newton iteration from the solution algorithm, and consequently reduces computational time. The algorithm presented in Table 2 was implemented into

ABAQUS / STANDARD

[40]

within the user-defined subroutine UMAT. Incompressibility was enforced in the FE simulation through a penalty method, where 1/D1 in (31) serves as the analogue of the bulk modulus. To maintain the incompressibility constraint, D1 was set to 10−9 Pa−1 for all simulations.

18

ACCEPTED MANUSCRIPT

Table 2: Algorithmic flowchart for damage–induced isotropic growth implemented in ABAQUS / STANDARD within the user-defined subroutine UMAT, [40].

ˆ iso (1) Given: F, Ψ max (2) Check growth criterion: D − Dth > 0  (3) Calculate growth stretch: λg = exp (k/3α2 )[1 − (1 + αt) exp(−αt)]hD − Dth i

(4) Calculate elastic deformation gradient: Fe = F/λg

CR IP T

(5) Calculate elastic right Cauchy–Green tensor: Ce = FT e Fe

(6) Calculate elastic second Piola–Kirchhoff stress tensor: Se = 2 (7) Calculate second Piola–Kirchhoff stress tensor: S = Se /λ2g

ˆ ∂Ψ ∂Ce

AN US

(8) Push–forward to Kirchhoff stress tensor: τ = FSFT ˆ ∂Se ∂ 2Ψ (9) Calculate elastic elasticity tensor: Ce = 2 =4 ∂Ce ∂Ce 2 (10) Push–forward to Eulerian elasticity tensor: (C)ijkl = Fe iI Fe jJ Fe kK Fe lL (Ce )IJKL

(11) Calculate co-rotational elasticity tensor: J(Λ)ijkl = J(C)ijkl +(τil δjk + τjk δil + τik δjl + τjl δik )/2

M

3.4. Verification

To verify the implemented model, we simulate the growth of a block of biological tissue sub-

ED

ject to uniaxial extension, and compare the FE results with the corresponding (semi-analytical) solutions obtained from MATLAB [52]. For all simulations we use a single 8-node linear brick,

PT

hybrid, constant pressure element (C3D8H). This hybrid formulation is used to prevent locking due to incompressibility. The material parameters are taken from Table 1. The block is first

CE

extended (loaded) up to a stretch of 1.4, and afterwards the block is unloaded to 1.0, assuming that damage occured during that loading path. Subsequently, the block grows over time due to induced damage. In Fig. 2 the evolution of the growth stretch λg is shown indicating perfect

AC

agreement between the finite element and the MATLAB results.

4. Parameter study of the growth model In this section, the sensitivity of the proposed model for damage-induced growth in arterial tissues is investigated. In particular we investigate the damage parameters r1 and m1 , and the 19

ACCEPTED MANUSCRIPT

1.05 1.045

λmax = 1.044 g

1.035 Matlab result

1.03

Finite element result

1.025 1.02 1.015 1.01 1.005 1 0

10

20

30

40

50

60

Time (days)

70

80

CR IP T

Growth stretch λg

1.04

90

100

Figure 2: Evolution of the growth stretch λg comparing finite element results with (semi-analytical) results

AN US

obtained from MATLAB indicating perfect agreement. The block is first stretched up to 1.4, and afterwards unloaded to 1.0 assuming that damage occured during that loading path. The material parameters are taken from Table 1.

growth parameters α and k in more detail. By considering growth of a block of biological

M

tissue subject to uniaxial extension, as described in the previous section, and by employing the parameters from Table 1, each time one of the damage parameters r1 , m1 and the growth

ED

parameters α, k are varied while the other parameters are kept fixed. Figure 3 shows plots between the growth stretch λg and the time up to 50 days, where the solid curves are produced with the parameters from Table 1, and the broken curves are obtained by multiplying the

PT

individual parameter with a factor ranging between 1.5 and 2.0 for r1 and 0.5 and 1.2 for m1 , α and k. The results indicate that the growth stretch λg is less sensitive to the variation of the

CE

parameter m1 amongst the four investigated parameters. For example, by increasing r1 , m1 , α and k by a factor of 1.5 (+50%), it results in −1.72%, −0.48%, −2.39% and 2.20% variation

AC

in the maximum growth stretch λg , respectively. Figures 3(c), (d) show that a reduction of α

entails an increase in the maximum growth stretch λg , while a reduction of k causes a decrease in the maximum λg . In addition, by increasing the parameters α and k the growth stretch λg reaches its maximum earlier and later, respectively.

20

ACCEPTED MANUSCRIPT

Growth stretch λg

1.04

1.05

r1 1.2r1 1.5r1 1.8r1 2.0r1

m1 1.2m1

Growth stretch λg

1.05

1.03

1.02

1.04

1.5m1 0.8m1

1.03

0.5m1

1.02

1.01

1.01

5

10

15

20

25

30

35

40

45

50

1 0

5

10

Time (days)

1.08

α 1.2α

k 1.07

Growth stretch λg

0.8α 0.5α

1.1

1.05

1.06 1.05

1.2k 1.5k 0.8k 0.5k

1.04

AN US

Growth stretch λg

25

30

35

40

45

50

(b)

1.5α 1.15

20

Time (days)

(a) 1.2

15

CR IP T

1 0

1.03 1.02 1.01

1 0

5

10

15

20

25

30

35

40

Time (days)

45

50

1 0

5

10

15

20

25

30

35

40

45

50

Time (days)

(c)

(d)

M

Figure 3: Parameter study of damage-induced growth of a block of biological tissue subject to uniaxial extension, as described in Section 3.4, showing relations between the growth stretch λg and the time up to 50 days. The

ED

damage parameters r1 and m1 and the growth parameters α and k are taken from Table 1; related results are shown by solid curves. Each individual parameter is varied by a multiplication factor, while the other parameters

PT

are kept fixed; related results are shown by broken curves.

CE

5. Finite element simulation of restenosis after angioplasty This example deals with the simulation of restenosis of an arterial wall occurring after angioplasty triggered by damage. The particular aim of this example is to show the general features

AC

of the proposed damage-induced growth model, and that the algorithm, as outlined in Table 2, works within the context of the FE method. This example is of an academic nature, and hence does not reflect a patient-specific geometry and measured material properties. The 3D FE model was built in

ABAQUS / STANDARD

where the initial geometry of the artery

is considered as a three-layered cylindrical segment with an inner diameter of 4 mm, a wall 21

ACCEPTED MANUSCRIPT

Table 3: Values for the hyperelastic, damage and growth parameters for the individual arterial layers (intima, media, adventitia) used for the simulation of restenosis after angioplasty.

Media

Adventitia

µ [MPa]

0.049

0.020

0.016

k1 [MPa]

15.467

0.180

0.845

k2 [–]

2.085

100

22.300

β [◦ ]

43.9

5.76

56.3

κ [–]

0.23

0.314

0.32

r1 [–]

1.37

3.36

3.29

0.0191

0.045

m1 [MPa] 0.0198

CR IP T

Intima

0.1

0.01

0.01

α [day−1 ]

0.05

0.05

0.05

k [day−2 ]

0.005

0.01

0.01

AN US

Dth [–]

thickness of 0.75 mm, and a length of 60 mm. The thickness ratios adventitia : media : intima

M

are assumed to be 2 : 3 : 1, and because of the underlying symmetry only one quarter of the arterial segment was discretized. In total 6 elements were used through-the-thickness of the

ED

arterial wall (2, 3, and 1 for adventitia, media, and intima, respectively). Eight-node linear brick, hybrid, constant pressure elements (C3D8H) were assigned to the mesh. We applied

PT

symmetry boundary conditions on two surfaces. One end of the quarter of the artery was free while at the other end the axial displacements were restricted, see Fig. 4.

CE

To analyze the problem, two steps were considered in

ABAQUS / STANDARD :

(i) ‘loading/

unloading’ during which damage occurred, and (ii) ‘growth’. To simulate balloon inflation/

AC

deflation a pressure of 40 kPa was applied at the inner surface of the artery (inflation) and then removed (deflation), i.e. the ‘loading/unloading’ step. The damage which was induced in the wall during that first step is the trigger for growth to build the neointima, i.e. within the second step of the simulation. The used values for the hyperelastic, damage and growth parameters for the individual arterial layers (intima, media, adventitia) are summarized in Table 3.

22

ACCEPTED MANUSCRIPT

CR IP T

Symmetry boundary condition

AN US

Axial constraint

Figure 4: Boundary conditions of a quarter of the artery for the finite element simulation of restenosis after angioplasty. ‘Symmetry boundary condition’ means that the circumferential displacements as well as the rotation about the axial and radial directions are restricted on both surfaces, and ‘axial constraint’ means that the displacements in the axial direction are fixed, while the other degrees of freedom are released.

M

Figure 5(a) depicts the FE results immediately after angioplasty, while Fig. 5(b) shows the results after the procedure at 180 days. The results show an increase in the thickness of the

ED

intimal layer due to the formation of a neointima. The evolution of the growth stretch λg with time (in days), for points on the inner surface of the intimal layer, is presented in Fig. 6. The

PT

results indicate that the growth is almost saturated after about 125 days and no significant growth occurs afterwards. Furthermore, from Fig. 6 it can be seen that the rate λ˙ g of growth

CE

is increased during the first 40 days after angioplasty, and afterwards it is decreased. These general observations are in a qualitative agreement with clinical observations (see, e.g., [47]).

AC

In this example we have focussed on a simplified simulation, and several clinically relevant issues were ignored in order to keep the computational effort relatively low. Computational modeling of the angioplasty process with more realistic features including patient–specific geometries and properties (as in, e.g., [15]), inelastic mechanical properties of the arterial wall including residual stresses (as in, e.g., [38]), and the effect of atherosclerotic plaques on the biomechanical response of the wall (as in, e.g., [53]) should be taken into consideration 23

CR IP T

ACCEPTED MANUSCRIPT

(a)

(b)

Figure 5: Finite element results of restenosis after angioplasty: (a) immediately after the procedure (after loading/unloading); (b) at 180 days after the procedure. Results show an increase in the thickness of the intimal layer

AN US

due to the formation of a neointima. 1.45 1.4

1.3 1.25

M

Growth stretch λg

1.35

1.2 1.15

1.05

20

PT

1 0

ED

1.1

40

60

80

100

120

140

160

180

Time (days)

(a)

CE

Figure 6: Evolution of the growth stretch λg over time (in days) after angioplasty due to damage for points on the inner surface of the intimal layer.

AC

in a more advanced study. However, the results presented here show the model capability to reproduce well-established features of in-stent restenosis.

6. Conclusion In the present study we proposed a microstructurally–based damage–induced growth model 24

ACCEPTED MANUSCRIPT

which we used to simulate in-stent restenosis. The finite growth theory was adopted and the multiplicative decomposition of the deformation gradient into an elastic and a growth part employed. The model was developed based on the changes that cells and tissues undergo after stent implementation, including vascular smooth muscle cell migration and proliferation.

happens on the cellular level after stent deployment.

CR IP T

The mechanobiological growth model provides a link between in-stent restenosis and what

By implementing the growth model into the nonlinear FE software

ABAQUS / STANDARD

within the user-defined subroutine UMAT, we analyzed growth of a block of biological tissue subject to uniaxial extension, and the implemented model has been verified by comparing the FE results with the corresponding (semi-analytical) solutions obtained from MATLAB. A pa-

AN US

rameter study of the model has also been conducted, which gives a clear picture of the effect of the involved parameters on the amount of growth.

Finally, the problem of restenosis after angioplasty was analyzed on the basis of a simplified example. It has been shown that the results obtained from the growth model are in qualitative

M

agreement with clinical observations. It should be pointed out that the simulation of restenosis after angioplasty, which was studied here, has a more or less academic character, because the

ED

used material parameters, the geometry and the boundary conditions were assumed. However, several features such as nonlinearity, anisotropy and the multi-layered structure of the artery wall were included. Several issues remain, and more complexity need to be considered to ob-

PT

tain more realistic and clinically relevant results. Nevertheless, in the presented FE example,

CE

the performance of the model was examined.

AC

Appendix A. Derivation of eq. (24) We start by writing eq. (3) in the index notation as −1 Ce ij = Fg−1 ni Cnm Fg mj .

Differentiation of (A.1) with respect to Ckl gives   ∂Ce ∂Ce ij ∂Cnm −1 −1 −1 −1 = = Fg−1 Fg mj = Fg−1 ni ni δnk δml Fg mj = Fg ki Fg lj . ∂C ijkl ∂Ckl ∂Ckl 25

(A.1)

(A.2)

ACCEPTED MANUSCRIPT

(A.3)

CR IP T

On use of (A.2)4 , the first term of eq. (23) can be represented as       ∂Ψ ∂Ce ˙ ∂Ψ ∂Ce ∂Ψ ˙ : :C= Ckl = F −1 F −1 C˙ kl ∂Ce ∂C ∂Ce ij ∂C ijkl ∂Ce ij g ki g lj     ∂Ψ −1 −1 ˙ −1 ∂Ψ −T = Fg ki F Ckl = Fg F C˙ kl ∂Ce ij g lj ∂Ce g kl   −1 ∂Ψ −T ˙ F : C. = Fg ∂Ce g

Now, we elaborate on the second term of eq. (23). Differentiation of (A.1) with respect to Fg kl gives ∂Ce ∂Fg



= ijkl

∂Fg−1 ∂Fg−1 ∂Ce ij ni mj −1 = Cnm Fg−1 . + F C mj g ni nm ∂Fg kl ∂Fg kl ∂Fg kl

AN US



(A.4)

It is noted that in eq. (23), C and Fg are considered as two independent variables so that ∂Cnm /∂Fg kl = 0. Then, the two partial derivatives in (A.4) are obtained by using the following identity [49] ∂A−1 ∂A



ijkl


∂A−1 1 −1 −1 ij −1 =− = Aik Alj + A−1 A il kj . ∂Akl 2

(A.5)

M



(A.6)

PT

ED

Thus, (A.4)2 can be represented as   1 ∂Ce −1 −1 −1 −1 −1 = − Ce zu (Fg−1 nk Fg li Fg zn Fg um Fg mj + Fg nl Fg ki Fg zn Fg um Fg mj ∂Fg ijkl 2 −1 −1 −1 −1 −1 +Fg−1 ni Fg mk Fg lj Fg zn Fg um + Fg ni Fg ml Fg kj Fg zn Fg um ),

AC

CE

where Cnm = Fg zn Ce zu Fg um has been used. After simplification of (A.6) we get  
1 ∂Ce −1 −1 −1 + δ δ F + δ δ F + δ δ F = − Ce zu δzk δju Fg−1 zi ul zl ju zi ku li g ki g lj g kj ∂Fg ijkl 2
1 −1 −1 −1 = − Ce kj Fg−1 li + Ce lj Fg ki + Ce ik Fg lj + Ce il Fg kj . 2

26

(A.7)

ACCEPTED MANUSCRIPT

CR IP T

Thus, the second term of eq. (23) can be written as    
∂Ψ ∂Ce ˙ ∂Ψ ∂Ce : : Fg = F˙ g kl ∂Ce ∂Fg ∂Ce ij ∂Fg ijkl "    ∂Ψ 1 ∂Ψ −1 ˙ ˙ Ce kj Fg li Fg kl + Ce lj Fg−1 =− ki Fg kl 2 ∂Ce ij ∂Ce ij #     ∂Ψ ∂Ψ ˙ ˙ + Ce ik Fg−1 Ce il Fg−1 lj Fg kl + kj Fg kl ∂Ce ij ∂Ce ij      
∂Ψ 1 ∂Ψ T −1 ˙ Ce =− Fg Fg ki + Ce F˙ g F−1 g 2 ∂Ce ki ∂Ce li li #      
∂Ψ ∂Ψ T + Ce F˙ g F−1 + Ce F˙ g F−1 , g g kj ∂Ce kj ∂Ce lj lj

(A.8)

where we have used the fact that ∂Ψ/∂Ce and Ce are symmetric tensors. We note that the

AN US

rotation is lumped into the elastic part of the deformation gradient and Fg is symmetric (see the Remark after eq. (46)). With this in mind, (A.8) can finally be represented as ∂Ψ ∂Ce ˙ ∂Ψ ˙ −1 : : Fg = −2Ce : Fg Fg . ∂Ce ∂Fg ∂Ce

(A.9)

With (10)1 and the property ∂Jg /∂Fg = Jg F−T g the third term in eq. (23) is obtained as

(A.10)

PT

ED

M

∂Ψ ∂ρ0g ˙ ∂Ψ ∂(ρg Jg ) ˙ ∂Ψ : Fg = : Fg = 0 ρg Jg F−T : F˙ g g 0 0 ∂ρg ∂Fg ∂ρg ∂Fg ∂ρg ∂Ψ 0 ∂Ψ I : F˙ g F−1 = ρ0g 0 tr(F˙ g F−1 g ) = ρg g . ∂ρg ∂ρ0g

CE

Appendix B. Derivation of the elasticity tensor (53)3 We start by writing (53)2 in the index notation as −1 ∂(Fg−1 ∂Se mn −1 im Se mn Fg jn ) = 2Fg−1 F im ∂Ckl g jn  ∂Ckl  ∂Se mn ∂Ce rq −1 ∂Ce rq −1 = Fg−1 Fg jn = Fg−1 F . im 2 im (Ce )mnrq ∂Ce rq ∂Ckl ∂Ckl g jn

AC

(C)ijkl = 2

(B.1)

In view of (A.2), eq. (B.1)4 can be reformulated as

−1 −1 −1 (C)ijkl = Fg−1 im Fg jn (Ce )mnrq Fg kr Fg lq .

Finally, (B.2) can be represented in the tensor notation as displayed in (53)3 . 27

(B.2)

ACCEPTED MANUSCRIPT

Appendix C. Derivation of the explicit expressions of the elasticity tensor This appendix deals with the calculation of the explicit expression for the isochoric part Ciso e m

of the elasticity tensor Ce . We start from (56) and consider first Ce which results from the matrix material. With the definition (57)1 and with (50)1 we obtain after some straightforward

m Ce

CR IP T

manipulations (see also [49]) ˜m) P:S e =2 ∂Ce ˜ m )P ˜ m : PT + 2 Tr(J −2/3 S ˜ − 2 (C−1 ⊗ Sm + Sm ⊗ C−1 ), =P:C e e e e e e 3 3 e −2/3

∂(Je

(C.1)

˜ = C−1 C−1 − C−1 ⊗ C−1 /3 is the where Tr(•) = (•) : Ce , is the trace of (•) and P e e e e

AN US

−4/3 ˜ m ˜m fourth-order projection tensor in the material description, while C ∂ Se /∂Ce is the e = 2Je

(fictitious) elasticity tensor in the material description. Next, we note that according to (50)2 , ˜ m is a zero tensor, and (C.1) can be simplified as C e

4 m ˜ − 2 (C−1 ⊗ Sm + Sm ⊗ C−1 ), Ce = Je−2/3 µI1e P e e e 3 3 e m

ED

M

where I1e = trCe . In view of (50), Se is obtained as     1 −1 1 m −1 −2/3 −2/3 Se = 2µJe I − Ce ⊗ Ce : I = 2µJe . I − I1e Ce 3 3

(C.2)

(C.3)

f,i

Second we calculate Ce , i.e. (57)2 , which relates to the two families of collagen fibers. We

  ˜ f,i = 2k1 E¯i exp(k2 E¯ 2 ) κI + (1 − 3κ)Ai , S e i 0

AC

and

CE

PT

˜ f,i and Sf,i can be obtained as note that, in view of (51) and (52), the explicit expression for S e e (C.4)

1 f,i i Se = 2Je−2/3 k1 E¯i exp(k2 E¯i2 )(I − C−1 e ⊗ Ce ) : [κI + (1 − 3κ)A0 ] 3   1 i −1 −2/3 2 ¯ ¯ = 2Je k1 Ei exp(k2 Ei ) κI + (1 − 3κ)A0 − [κI1e + (1 − 3κ)Iie ]Ce , (C.5) 3

in which Iie = ai0 · Ce ai0 , i = 4, 6, has been utilized.

28

ACCEPTED MANUSCRIPT

f,i

Then, by employing the same procedure which let to (C.1)2 , from (57)2 we may deduce Ce by using (C.4) and (C.5)2 . Thus,

CR IP T

2 2 f,i f,i f,i −1 T −2/3 ˜ f,i ˜ ˜ f,i Ce = P : C Se )P − (C−1 e : P + Tr(Je e ⊗ Se + Se ⊗ Ce ) 3 3 4 −2/3 ¯ f,i T ˜ ˜ = P : Ce : P + Je k1 Ei exp(k2 E¯i2 )[κI1e + (1 − 3κ)Iie ]P 3  4 −1 − Je−2/3 k1 E¯i exp(k2 E¯i2 ) κ(C−1 e ⊗ I + I ⊗ Ce ) 3 2 i i −1 −1 −1 +(1 − 3κ)(C−1 e ⊗ A0 + A0 ⊗ Ce ) − [κI1e + (1 − 3κ)Iie ]Ce ⊗ Ce }, (C.6) 3 ˜ f,i can be calculated using (52)2 , i.e. where, in view of (C.4), C e

AN US

˜ f,i ˜ f,i ¯ −4/3 ∂ Se f,i −4/3 ∂ Se ∂ Ei ˜ Ce = 2Je = 2Je ∂ E¯i ∂Ce ∂Ce = 4Je−4/3 k1 exp(k2 E¯i2 )(1 + 2k2 E¯i2 )

×[κ2 I ⊗ I + (1 − 3κ)2 Ai0 ⊗ Ai0 + κ(1 − 3κ)(I ⊗ Ai0 + Ai0 ⊗ I)], an expression which is needed in (C.6)2 .

M

References

(C.7)

downloaded from:

ED

[1] World Health Organization, Global status report on noncommunicable diseases 2014,

PT

http://www.who.int/nmh/publications/ncd-status-report-2014/en/

[2] Writing Group Members, D. Mozaffarian, E. J. Benjamin, A. S. Go, D. K. Arnett, M. J.

CE

Blaha, M. Cushman, S. T. Das, S. de Ferranti, J. P. Despr´es, H. J. Fullerton, V. J. Howard, M. D. Huffman, C. R. Isasi, M. C. Jim´enez, S. E. Judd, B. M. Kissela, J. H. Lichtman, L. D. Lisabeth, S. Liu, R. H. Mackey, D. J. Magid, D. K. McGuire, E. R. Mohler 3rd,

AC

C. S. Moy, P. Muntner, M. E. Mussolino, K. Nasir, R. W. Neumar, G. Nichol, L. Palaniappan, D. K. Pandey, M. J. Reeves, C. J. Rodriguez, W. Rosamond, P. D. Sorlie, J. Stein, A. Towfighi, T. N. Turan, S. S. Virani, D. Woo, R. W. Yeh, M. B. Turner, American Heart Association Statistics Committee, Stroke Statistics Subcommittee, Executive summary: Heart disease and stroke statistics – 2016 update: A report from the american heart association, Circulation 133 (2016) 447–454. 29

ACCEPTED MANUSCRIPT

[3] G. A. Holzapfel, T. C. Gasser, Computational stress–deformation analysis of arterial walls including high-pressure response, Int. J. Cardiol. 116 (2007) 78–85. [4] C. J. Boyle, A. B. Lennon, P. J. Prendergast, In silico prediction of the mechanobiological response of arterial tissue: Application to angioplasty and stenting, J. Biomech. Eng. 133

CR IP T

(2011) 081001. [5] G. Dangas, F. Kuepper, Restenosis: Repeat narrowing of a coronary artery prevention and treatment, Circulation 105 (2002) 2586–2587.

[6] R. Hoffmann, G. Mintz, Coronary in-stent restenosis – predictors, treatment and preven-

AN US

tion, Eur. Heart J. 21 (2000) 1739–1749.

[7] J. Daemen, P. W. Serruys, Drug-eluting stent update 2007: part II: Unsettled issues, Circulation 116 (2007) 961–968.

[8] H. Tahir, A. G. Hoekstra, E. Lorenz, P. V. Lawford, D. R. Hose, J. Gunn, D. J. W.

M

Evans, Multi-scale simulations of the dynamics of in-stent restenosis: Impact of stent deployment and design, Interface Focus 1 (2011) 365–373.

ED

[9] H. Tahir, I. Niculescu, C. Bona-Casas, R. M. Merks, A. G. Hoekstra, An in silico study on the role of smooth muscle cell migration in neointimal formation after coronary stent-

PT

ing, J. R. Soc. Interface 12 (2015) 20150358. [10] R. Kornowski, M. K. Hong, F. O. Tio, O. Bramwell, H. Wu, M. B. Leon, In-stent resteno-

CE

sis: contributions of inflammatory responses and arterial injury to neointimal hyperplasia, J. Am. Coll. Cardiol. 31 (1998) 224–230.

AC

[11] H. Wieneke, M. Haude, M. Knocks, A. Gutersohn, C. Von Birgelen, D. Baumgart, R. Erbel, Evaluation of coronary stents in the animal model: A review, Materialwissenschaft und Werkstofftechnik 30 (1999) 809–813.

[12] F. G. Welt, C. Rogers, Inflammation and restenosis in the stent era, Arterioscler. Thromb. Vasc. Biol. 22 (2002) 1769–1776. 30

ACCEPTED MANUSCRIPT

[13] A. K. Mitra, D. K. Agrawal, In stent restenosis: Bane of the stent era, J. Clin. Pathol. 59 (2006) 232–239. [14] F. Migliavacca, L. Petrini, M. Colombo, F. Auricchio, R. Pietrabissa, Mechanical behavior of coronary stents investigated through the finite element method, J. Biomech. 35

CR IP T

(2002) 803–811. [15] G. A. Holzapfel, M. Stadler, T. C. Gasser, Changes in the mechanical environment of stenotic arteries during interaction with stents: Computational assessment of parametric stent design, J. Biomech. Eng. 127 (2005) 166–180.

[16] C. Lally, F. Dolan, P. J. Prendergast, Cardiovascular stent design and vessel stresses: a

AN US

finite element analysis, J. Biomech. 38 (2005) 1574–1581.

[17] W. Q. Wang, D. K. Liang, D. Z. Yang, M. Qi, Analysis of the transient expansion behavior and design optimization of coronary stents by finite element method, J. Biomech. 39 (2006) 21–32.

M

[18] M. De Beule, P. Mortier, S. G. Carlier, B. Verhegghe, R. Van Impe, P. Verdonck, Re-

(2008) 383–389.

ED

alistic finite element-based stent design: the impact of balloon folding, J. Biomech. 41

PT

[19] P. Mortier, G. A. Holzapfel, M. De Beule, D. Van Loo, Y. Taeymans, P. Segers, P. Verdonck, B. Verhegghe, A novel simulation strategy for stent insertion and deployment

CE

in curved coronary bifurcations: comparison of three drug-eluting stents, Ann. Biomed. Eng. 38 (2010) 88–99.

AC

[20] H. Zahedmanesh, D. J. Kelly, C. Lally, Simulation of a balloon expandable stent in a realistic coronary artery – determination of the optimum modelling strategy, J. Biomech. 43 (2010) 2126–2132.

[21] L. H. Timmins, M. R. Moreno, C. A. Meyer, J. C. Criscione, A. Rachev, J. E. Moore Jr., Stented artery biomechanics and device design optimization, Med. Biol. Eng. Comput. 45 (2007) 505–513. 31

ACCEPTED MANUSCRIPT

[22] H. Zahedmanesh, C. Lally, Determination of the influence of stent strut thickness using the finite element method: implications for vascular injury and in-stent restenosis, Med. Biol. Eng. Comput. 47 (2009) 385–393. [23] I. Pericevic, C. Lally, D. Toner, D. J. Kelly, The influence of plaque composition on underlying arterial wall stress during stent expansion: the case for lesion-specific stents,

CR IP T

Med. Eng. Phys. 31 (2009) 428–433.

[24] G. A. Holzapfel, J. J. Mulvihill, E. M. Cunnane, M. T. Walsh, Computational approaches for analyzing the mechanics of atherosclerotic plaques: a review, J. Biomech. 47 (2014) 859–869.

AN US

[25] F. Migliavacca, L. Petrini, P. Massarotti, S. Schievano, F. Auricchio, G. Dubini, Stainless and shape memory alloy coronary stents: a computational study on the interaction with the vascular wall, Biomech. Model. Mechanobiol. 2 (2004) 205–217. [26] H. Zahedmanesh, C. Lally, A multiscale mechanobiological model using agent based

M

models; application to vascular tissue engineering, Biomech. Model. Mechanobiol. 11

ED

(2012) 363–377.

[27] A. Masselot, B. Chopard, A lattice Boltzmann model for particle transport and deposi-

PT

tion, EPL (Europhysics Letters) 42 (1998) 259–264. [28] A. Ilachinski, Cellular Automata: A Discrete Universe, World Scientific Publishing Co.,

CE

River Edge, NJ, 2001.

[29] D. C. Walker, J. S. Southgate, G. Hill, M. Holcombe, D. R. Hose, S. M. Wood,

AC

S. Mac Neil, R. H. Smallwood, The epitheliome: agent-based modelling of the social behaviour of cells, Biosystems 76 (2004) 89–100.

[30] M. Wooldridge, An Introduction to MultiAgent Systems, 2nd Edition, John Wiley & Sons, 2009.

32

ACCEPTED MANUSCRIPT

[31] D. J. Evans, P. V. Lawford, J. Gunn, D. Walker, D. R. Hose, R. H. Smallwood, B. Chopard, M. Krafczyk, J. Bernsdorf, A. Hoekstra, The application of multiscale modelling to the process of development and prevention of stenosis in a stented coronary artery, Philos. Trans. A Math. Phys. Eng. Sci. 366 (2008) 3343–3360. [32] E. Adiguzel, P. J. Ahmad, C. Franco, M. P. Bendeck, Collagens in the progression and

CR IP T

complications of atherosclerosis, Vasc. Med. 14 (2009) 73–89.

[33] F. J. Gijsen, F. Migliavacca, S. Schievano, L. Socci, L. Petrini, A. Thury, J. J. Wentzel, A. F. van der Steen, P. W. Serruys, G. Dubini, Simulation of stent deployment in a realistic human coronary artery, Biomed. Eng. Online 7 (2008) 23.

AN US

[34] G. Himpel, Computational modeling of biomechanical phenomena: Remodeling, growth and reorientation, PhD dissertation, University of Kaiserslautern, Germany (2007). [35] S. G¨oktepe, O. J. Abilez, E. Kuhl, A generic approach towards finite growth with ex-

Solids 58 (2010) 1661–1680.

M

amples of athlete’s heart, cardiac dilation, and cardiac wall thickening, J. Mech. Phys.

ED

[36] E. K. Rodriguez, A. Hoger, A. D. McCulloch, Stress-dependent finite growth in soft elastic tissues, J. Biomech. 27 (1994) 455–467.

PT

[37] V. A. Lubarda, A. Hoger, On the mechanics of solids with a growing mass, Int. J. Solids Structures 39 (2002) 4627–4664.

CE

[38] B. Fereidoonnezhad, R. Naghdabadi, G. A. Holzapfel, Stress softening and permanent deformation in human aortas: continuum and computational modeling with application

AC

to arterial clamping, J. Mech. Behav. Biomed. Mater. 61 (2016) 600–616.

[39] T. C. Gasser, R. W. Ogden, G. A. Holzapfel, Hyperelastic modelling of arterial layers with distributed collagen fibre orientations, J. R. Soc. Interface 3 (2006) 15–35.

[40] Abaqus 6.13-4 Analysis User’s Guide, Dassault Syst`emes Simulia Corp., 2013.

33

ACCEPTED MANUSCRIPT

[41] FEAP – A Finite Element Analysis Program, Version 8.2 User Manual, University of California at Berkeley, Berkeley, California, 2008. [42] J. Helfenstein, M. Jabareen, E. Mazza, S. Govindjee, On non-physical response in models for fiber-reinforced hyperelastic materials, Int. J. Solids Struct., 47 (2010) 2056–

CR IP T

2061. [43] W. Ehlers, G. Eipper, The simple tension problem at large volumetric strains computed from finite hyperelastic material laws, Acta Mech. 130 (1998) 17–27.

[44] G. A. Holzapfel, J. A. Niestrawska, R. W. Ogden, A. J. Reinisch, A. J. Schriefl, Modelling non-symmetric collagen fibre dispersion in arterial walls, J. R. Soc. Interface 12

AN US

(2015) 20150188.

[45] J. A. Niestrawska, Ch. Viertler, P. Regitnig, T. U. Cohnert, G. Sommer, G. A. Holzapfel, Microstructure and mechanics of healthy and aneurysmatic abdominal aortas: experi-

M

mental analysis and modeling, J. R. Soc. Interface, in press.

[46] R. Fukui, M. Amakawa, M. Hoshiga, N. Shibata, E. Kohbayashi, M. Seto, Y. Sasaki,

ED

T. Ueno, N. Negoro, T. Nakakoji, M. Ii, F. Nishiguchi, T. Ishihara, N. Ohsawa, Increased migration in late G(1) phase in cultured smooth muscle cells, Am. J. Physiol.

PT

Cell Physiol. 279 (2000) C999–1007. [47] R. S. Schwartz, A. Chu, W. D. Edwards, S. S. Srivatsa, R. D. Simari, J. M. Isner, D. R.

CE

Holmes Jr., A proliferation analysis of arterial neointimal hyperplasia: lessons for antiproliferative restenosis therapies, Int. J. Cardiol. 53 (1996) 71–80.

AC

[48] H. Weisbecker, D. M. Pierce, P. Regitnig, G. A. Holzapfel, Layer-specific damage experiments and modeling of human thoracic and abdominal aortas with non-atherosclerotic intimal thickening, J. Mech. Behav. Biomed. Mater. 12 (2012) 93–106.

[49] G. A. Holzapfel, Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons, Chichester, 2000.

34

ACCEPTED MANUSCRIPT

[50] A. Javili, B. Dortdivanlioglu, E. Kuhl, C. Linder, Computational aspects of growthinduced instabilities through eigenvalue analysis, Comput. Mech. 56 (2015) 405–420. [51] J. M. Young, J. Yao, A. Ramasubramanian, L. A. Taber, R. Perucchio, Automatic generation of user material subroutines for biomechanical growth analysis, J. Biomech. Eng.

CR IP T

132 (2010) 104505. [52] MATLAB, R2016a, The MathWorks Inc., Natick, MA, USA, 2016.

[53] D. E. Kiousis, S. F. Rubinigg, M. Auer, G. A. Holzapfel, A methodology to analyze changes in lipid core and calcification onto fibrous cap vulnerability: The human atherosclerotic carotid bifurcation as an illustratory example, J. Biomech. Eng. 131

AC

CE

PT

ED

M

AN US

(2009) 121002.

35