PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Ann Biomed Eng. Author manuscript; available in PMC 2016 July 1.
Published in final edited form as:
PMCID: PMC4497867
NIHMSID: NIHMS669011

Coupled simulation of hemodynamics and vascular growth and remodeling in a subject-specific geometry

Abstract

A computational framework to couple vascular growth and remodeling (G&R) with blood flow simulation in a 3D patient-specific geometry is presented. Hyperelastic and anisotropic properties are considered for the vessel wall material and a constrained mixture model is used to represent multiple constituents in the vessel wall, which was modeled as a membrane. The coupled simulation is divided into two time scales–a longer time scale for G&R and a shorter time scale for fluid dynamics simulation. G&R is simulated to evolve the boundary of the fluid domain, and fluid simulation is in turn used to generate wall shear stress and transmural pressure data that regulates G&R. To minimize required computation cost, the fluid dynamics are only simulated when G&R causes significant vascular geometric change. For demonstration, this coupled model was used to study the influence of stress-mediated growth parameters, and blood flow mechanics, on the behavior of the vascular tissue growth in a model of the infrarenal aorta derived from medical image data.

Keywords: aneurysm, blood flow, constrained mixture, coupled simulation, growth and remodeling

1 Introduction

Cardiovascular disease is one of the leading causes of morbidity and mortality in the world; the progression of which is closely related to vascular adaptation in response to mechanical and chemical stimuli. As a framework to better understand this adaptive behavior, a theory for vascular growth and remodeling (G&R) was proposed by Humphrey, et al., based on a constrained mixture model [11] in which the vessel wall is assumed to have the ability to adapt to vascular stimuli and recover a homeostatic state via smooth muscle cell synthesis and matrix turnover [22, 23]. Kinematic models of G&R have also been considered [19], where the time rates of change of the growth stretch ratios are assumed to linearly depend on the local smooth muscle stress and on vascular wall shear stress. The constrained mixture theory of G&R has been applied to study stress-mediated aneurysm expansion in idealized ellipsoidal and cylindrical geometry [2, 3] to better understand key factors influencing expansion rate and resulting aneurysm shape. Recently, G&R simulation has been extended to 3D geometries to predict the more complex case of asymmetric expansion [32]. In other works, [29] studied the role of collagen properties on AAA progression modeling, and [28] studied the influence of the initial state of the aorta on enlargement and mechanical behavior. Humphrey and Holzapfel [10] provided a review of the experimental data, computational models, mechanobiological factors, and open problems for G&R of human abdominal aortic aneurysm.

Vascular growth and remodeling is an inherently coupled problem involving both the dynamics of the vessel wall and stresses transmitted by blood flow. Recently, several papers have focused on extending G&R simulations to include coupling with blood flow mechanics. In [6], a theory of “small-on-large” was used to couple vascular G&R and hemodynamic simulation in a cylindrical geometry, where hemodynamic forces acting on the vessel wall were obtained from fluid-solid-interaction. Cylindrical sections have also been considered in coupled models of Watton, et al., [27]. Because both structural and fluid mechanics depend critically on vascular morphology, which is highly subject-specific, there has remained compelling need to extend coupled G&R simulations to fully 3D subject-specific geometries. Challenges in implementing such modeling include coupling simulation of the different time scales involved, defining anisotropic material properties in a subject-specific geometry, as well as the normal difficulties of simulating fluid flow and structural mechanics in complex vascular domains. An important step in this direction was implemented in [18], where G&R was applied to a cylindrical portion of an anatomically-realistic abdominal aorta, and simulation of the two time scales was obtained in a stepwise manner.

In this paper, we follow the basic framework of the constrained mixture theory for G&R described in [3, 21, 24]. Unlike previous works, we extend the formulation of coupling hemodynamics with vascular G&R to a realistic vascular model that includes branching and several arterial segments. This provides capability to consider more complete patient-specific geometries than previously considered, which is necessary to translate such modeling to many important clinical applications. The model used herein is derived from medical image data and vascular adaptation resulting from an idealized injury model is simulated. Blood flow is modeled as an incompressible, Newtonian fluid and simulated when G&R introduces significant geometric change. Hemodynamic forces are used to regulate G&R over longer time scales. In Sec. 2, we define the basic concepts, G&R kinetics, constitutive relations, hemodynamics and stress mediated growth laws. An algorithm for the coupled simulation is presented, which was implemented within a finite element framework using custom code and COMSOL as a generic FEM solver. In Sec. 3, we present simulations of different scenarios for abdominal aortic aneurysm expansion by varying initial mass loss. Coupled simulations are performed with different values of G&R feedback gains to demonstrate the influence of these constants on both the resulting aneurysm shape and hemodynamics.

2 Methods

2.1 Definitions

Vascular G&R is modeled by the theoretical framework proposed in [11]. The vessel wall is modeled as a membrane and treated as a constrained mixture, which implies that at each location the mixture (collagen + elastin + smooth muscle ) deforms together. The reference configuration will be denoted κ0 for the mixture. This configuration corresponds to zero transmural pressure, P = 0, however, constituents in this configuration are not necessarily stress-free due to prestress. At any time t, the deformed configuration of the mixture is defined as κ(t). The deformation gradient tensor F(t) maps κ0 [mapsto] κ(t).

Collagen is assumed to be anisotropic and is characterized by discretized collagen families k, which have individual collagen fiber directions and reference configurations. Collagen can be incorporated into the mixture at intermediate times τ, with pre-stretch defined by tensor Gk(τ). For each orientation k, Gk(τ) maps the natural configuration of the newly produced collagen at time τ, denoted κn(τ)k, to the deformed configuration of the mixture at this time τ; subscript n(τ) denotes the natural configuration at time τ. For each collagen family, we define a deformation gradient tensor, which maps the natural configuration of that family to the current deformed mixture configuration, by

Fn(τ)k(t)=F(t)F-1(τ)Gk(τ),
(1)

hence, the right Cauchy-Green deformation tensor can be obtained as

Cn(τ)k(t)=Fn(τ)k(t)Fn(τ)k(t).
(2)

These mappings and their compositions are shown schematically in Fig. 1.

Fig. 1
Configurations and associated mappings used to describe G&R framework.

The unit vector in the direction of collagen family k at time t is denoted ek(t), and collagen fiber produced at time τ is assumed to be deposited in the ek(τ) direction, i.e., in the direction of existing collagen. Therefore, the direction of the collagen fiber produced at time τ in its natural configuration, κn(τ)k, is described by the unit vector

en(τ)k=Gk(τ)-1ek(τ)Gk(τ)-1ek(τ).
(3)

With ek(τ) and en(τ)k, Gk(τ) can be defined as a two-point tensor

Gk(τ)=Gkek(τ)en(τ)k,
(4)

where Gk is the stretch ratio of newly produced collagen fiber from the natural configuration κn(τ)k to the mixture configuration κτ. We assume Gk does not depends on collagen family k and is equal to the homeostatic collagen stretch ratio Ghc, set herein to 1.05; subscript h denotes the homeostatic state. The stretch ratio Gk can be obtained as

λn(τ)k(t)=en(τ)k·Cn(τ)k(t)en(τ)k=F(t)F-1(τ)Gk(τ)en(τ)k=F(t)F-1(τ)Ghc(ek(τ)en(τ)k)en(τ)k=GhcF(t)F-1(τ)ek(τ)=GhcF(t)F-1(τ)ek(τ)/F-1(τ)ek(τ)ek(τ)/F-1(τ)ek(τ)=Ghcλk(t)λk(τ),
(5)

where λk(t) is the stretch ratio of the mixture in the direction of collagen family k, defined as

λk(t)=F(t)F-1(τ)ek(τ)F-1(τ)ek(τ).
(6)

For elastin, we assume the mapping from the natural configuration to the mixture reference configuration is Ge, and therefore the mapping from elastin’s natural configuration to the current mixture configuration is given by

Fne(t)=F(t)Ge.
(7)

In this model we take Ge=diag[G1e,G2e,1G1eG2e], which assumes the prestretch occurs in the principal directions of the existing mixture, where the third component has been written in terms of the prior two by imposing incompressibility.

2.2 Kinetics of growth & remodeling

The vessel wall has ability to adapt in response to mechanical stimuli in order to recover a homeostatic state. This process occurs by both the removal of old constituents, and incorporation of new constituents into the mixture due to natural turnover as well as stress-mediated growth. Let Mk(t) be the mass per unit area of each collagen family k, and the time evolution of Mk(t) is described by

Mk(t)=Mk(0)Qk(t)+0tmk(τ)qk(t-τ)dτ.
(8)

The first term on the right of the equation describes natural turnover of the initial mass produced before the G&R process. The second term describes the incorporation and natural turnover of the newly produced constituent. Qk(t) is the remaining fraction at time t for the kth collagen family produced at time 0. mk(τ) is the mass production rate of the kth collagen family at time τ and qk(tτ) is the remaining fraction at time t.

For elastin,

Me(t)=Me(0)Qe(t),
(9)

where Qe(t) is the corresponding remaining fraction of initial mass for elastin at time t. Since functional elastin is thought to be mainly produced during early development, there is no second term for newly produced elastin in (9), as there was for collagen in (8).

The remaining fraction for collagen family k is assumed to have the following form [6]

qk(t)={1,if0t<t112{cos(πt2-t1(t-t1))+1},ift1tt20,ift2<t,
(10)

where t2 is the lifespan of the constituent. For collagen we use t1/t2 = 0.2, and t2 = 1. Hence, time herein is normalized by the lifespan of collagen, which is generally between 70–80 days [29, 13]. The function representing the remaining fraction at time t of the initial mass is given by

Qk(t)=1-0tqk(τ)dτ0t2qk(τ)dτ=1-20tqk(τ)dτt1+t2.
(11)

2.3 Constitutive relations

Collagen was assumed to follow a Fung-type exponential constitutive relation. The strain energy per unit mass of the kth family of collagen is

Wk(In(τ)k)=c24c3{exp[c3(In(τ)k-1)2]-1},
(12)

where

In(τ)k=(λn(τ)k(t))2.
(13)

Let

Cn(τ)k(t)=Fn(τ)k(t)Fn(τ)k(t)
(14)

denote the right Cauchy-Green deformation tensor of the kth family of collagen produced at time τ, whose current direction is denoted by ek(t), with respect to its natural configuration. The subscript n(τ) denotes the natural configuration of the constituent produced at time τ. For elastin, the strain energy function is defined based on a Neo-Hookean behavior,

We(Cne(t))=c12(I1-3),
(15)

where

I1(t)=tr(Cne(t)).
(16)

Cne(t)=Fne(t)Fne(t) is the right Cauchy-Green deformation tensor of elastin with respect to its nature configuration.

Based on the mass-averaged principle for a constrained mixture, the total strain energy per unit area for all constituents at time t is

w(t)=we(t)+kwk(t)=Me(0)Qe(t)We(Cne(t))+k{Mk(0)Qk(t)Wk(In(0)k)+0tmk(τ)qk(t-τ)Wk(In(τ)k)}
(17)

where we(t) and wk(t) denote the total strain energy contributed by elastin and the kth collagen family, respectively. We do not directly model the passive response of smooth muscle due to the fact that smooth muscle is much more compliant than collagen and thus has minimal contribution to the passive mechanical behavior of the vessel wall [4]. However the active response of smooth muscle does play an important role in vascular adaptation to altered flow. It is important to note that the strain energy function depends on both the current state and the history of deformation. Once the mathematical form of the total strain energy function is formulated, the deformation of the vessel wall can be obtained from the virtual work principle

δI=SδwdA-sPn·δxda=0,
(18)

where P is the vascular transmural pressure, n is the normal vector on vessel wall surface and δx is the virtual displacement of the vessel wall. S and s denote the surface area of vessel wall in the reference and current configuration, respectively.

2.4 Defining local anisotropic material property

One of the challenges in implementing growth and remodeling in 3D patient specific geometry is to define local anisotropic material properties, i.e., to define the local directions for collagen families. Usually local collagen directions e0k are defined with respect to local circumferential and axial directions in the reference configuration, and later evolve with the mixture deformation F(t). Therefore, the current collagen directions ek(t) are defined as

ek(t)=F(t)e0k.
(19)

Note that normalization is needed to obtain unit direction.

In idealized geometries, circumferential and axial directions can be readily defined (see [2, 3]), but defining these directions in 3D patient specific geometry is nontrivial. In [31], the authors defined the local circumferential and axial directions based on a 2-D parameterization of the vessel wall surface. However, it is not obvious how to apply this approach to a geometry with multiple outlets or bifurcations. In the work herein, local circumferential and axial directions are defined using local principal curvature directions, from which, local collagen directions are then calculated. As long as the geometry is smooth enough (which is usually the case), principal curvature directions are everywhere well-defined.

2.5 Hemodynamics

Blood was modeled as an incompressible, Newtonian fluid described by the Navier-Stokes and continuity equation

ρb(v(x,t)t+v(x,t)·v(x,t))=-p(x,t)+μ2v(x,t),
(20)
·v(x,t)=0.
(21)

The blood density ρb and viscosity μ were set to 1.05 g/cm3 and 0.035 Poise. A no-slip, no-penetration boundary condition was specified along lumen surface of the fluid domain, and a velocity profile (Dirichlet boundary condition) was specified at the inlet plane. At the outflow surfaces, Neumann-type boundary conditions were specified by coupling resistive models of the downstream vascular beds. Namely, the pressure P0(t) at each outflow boundary was solved as

Po(t)=RQo(t)
(22)

The flow rate Qo at the respective outlets was obtained from the 3D domain, and using Eq. (22) the pressure Po at the outlet was computed and applied at the zero traction outlet faces as Neumman boundary conditions. Details of the FEM implementation of the boundary conditions can be found in [5].

It is important to note that the time scales for G&R ([less, similar]weeks) and blood flow simulation ([less, similar]second) are several orders of magnitude different. A fully-coupled simulation of both processes is neither efficient nor necessary, because the hemodynamics stresses imparted from the blood flow do not change significantly until usually several weeks of geometric change has occurred through G&R. Nominally, blood flow was simulated when G&R caused changes to the boundary of the fluid domain, and for all times in between, the values for WSS and pressure were well approximated by the values from the last blood flow simulation. Specifically, it was assumed that mesh quality was an appropriate measure to monitor vessel deformation and the need for computing a new Navier-Stokes solution. This is (at least) consistent with the fact that the numerical accuracy of the Navier-Stokes solution, and hence WSS, depends on mesh quality.

2.6 Stress-mediated growth & remodeling

Here we describe how vascular adaptation is regulated by wall tension and WSS. Once the deformation is obtained by solving the variational equation (18), the Cauchy stress tensor is obtained as

T(t)=1J(t)F(t)wI1(t)I1(t)F(t)+1J(t)F(t)kwIn(τ)k(t)In(τ)k(t)F(t)+Tactive=2J(t)weI1(t)B(t)+2J(t)kwkIn(τ)kIn(τ)kIk(t)ekek+Tactive,
(23)

where B = FF[top top] is the left Cauchy-Green deformation tensor, and J(t) is the determinant of the deformation gradient tensor F(t). The first term on the right hand side denotes the stress contribution of elastin and the second term denotes that of two collagen families. The last term Tactive denotes the active membrane stress due to active smooth muscle contraction and relaxation. However, to keep the model simple and computation tractable, we do not include this active response term in our later simulations. A scalar measure is obtained from the Cauchy stress tensor in the direction of collagen family k as

σk=ek·Tekh,
(24)

which is the stress in the ek direction. The thickness of the vessel wall was calculated as

h(t)=M(t)Jρ,
(25)

where M (t) = Σk Mk(t) is the total collagen mass and ρ denotes the volume density of collagen.

In G&R theory, the vascular homeostatic state is recovered through stress-mediated feedback. The mass production rate of the kth family of collagen is assumed to depend linearly on the deviation of wall tension σk with respect to the homeostatic value σh. In addition, WSS regulation may be active via interactions between vascular endothelium and blood flow. It is well-known that endothelial cells can release multiple vasoactive chemicals in response to altered WSS. These chemicals can affect smooth muscle cell proliferation and collagen turnover. For instance, nitric oxide (NO) is a potential inhibitor of synthesis of collagen and smooth muscle proliferation [15], while endothelin-1 (ET-1) is a promoter of synthesis of collagen and smooth muscle proliferation [14]. NO is up-regulated in response to increased WSS [20] and ET-1 is up-regulated by decreased WSS [12]. These factors make WSS regulation uncertain. Due to lack of further information, we naively assume mass production rate of collagen is proportional to the deviation of wall shear stress τw from the preferred state τwh and consider cases of both positive and negative feedback. Thus the complete stress mediated growth law is given by

mk(t)=M(t)M(0)(Kσ(σk(t)-σh)-Kτ(τw(t)-τwh)+fhk),
(26)

where fhk is the basal value of mass production rate for collagen family k,

fhk=Mk(0)0qk(τ)dτ,
(27)

which balances the degradation rate of collagen in the homeostatic state. Note that τw(t) in (26) is given directly by the blood flow simulation, whereas σk depends on transmural pressure. Kσ and Kτ are feedback gains for the deviations of wall tension and wall shear stress. For Kσ > 0, if σk(t) > σh, the mass production rate mk(t) will increase, and in turn, cause an increase in wall thickness h(t). The thickening of the vessel wall will decrease the vessel radius for the same transmural pressure. Based on Laplace’s law for a cylinder, σ=Prh, and we can postulate that σk will hence decrease and return back to the homeostatic value. This means that for Kσ > 0, the growth law (26) initiates a negative feedback mechanism. Using similar arguments for wall shear stress based on simple Poiseuille flow, τ=4ηQπr3, we can postulate that a decreased vessel radius will increase WSS, so that Kτ > 0 also initiates a negative feedback mechanism for wall shear stress deviation. Therefore, if Kσ and Kτ are large enough, the stress will converge to the corresponding homeostatic value.

3 Results

In this section we apply G&R simulation to a vascular model of the aorta whose lumen morphology is derived from medical image data. The model is a truncated section of an aortofemeral model (ID: OSMSC0006) publicly available on vascularmodel.org [30]. After initialization of the model, scenarios for G&R are considered by introducing mass loss at various locations to simulate the development of abdominal aortic aneurysm. Material constants used in the simulations are listed in Table 1. In the following results, two main helical collagen directions (45° and 135° from the axial direction) are included in the simulations. This is consistent with prior observations that close to 90% of the total mass of collagen is distributed in these helical directions [9], although this does not hold universally true and additional families or orientations could be considered. Fig. 3 displays the fiber orientation about the aortic bifurcation region. It is observed that the collagen directions in both families naturally become aligned about the apical ridge at the bifurcation, which is consistent with the simulation results in [8] and consistent with experimental findings in [7].

Fig. 3
Local collagen directions calculated based on principal curvature directions at the aortic bifurcation: (a) 45° collagen family, (b) 135° collagen family.
Table 1
Material constants for elastin and collagen [32, 25]

3.1 Generation of initial homeostatic configuration

In the G&R theory, the healthy vessel wall is assumed to be in the homeostatic state, i.e., the stress of each constituent is equal to the homeostatic value, and therefore no significant growth and remodeling is induced beyond the basal production and turnover rate. Recall that a membrane model is used for the vessel wall. The vascular geometry derived from the medical image data is irregular, and therefore achieving an initial homeostatic stress state requires specifying an appropriately varying mass distribution of the constituents. This can be accomplished by an initial G&R stage to solve for an appropriate homeostatic mass distribution, as proposed in [32]. Namely, an accelerated G&R stage is simulated using a feedback gain for wall tension, Kσ, set to 0.005. This enabled the mass distribution to converge (approximately) to the homeostatic state without significantly altering the overall geometry.

In Fig. 2(a)(b), and 2(c)(d), the deviation of stress relative to the homeostatic values for the 45°, and 135°, respectively, oriented collagen family are plotted before and after the initial homeostatic simulation. Wall tension distributions become mostly uniform after the initial homeostatic simulation, with stress values in both directions generally within the range of 0.9σh to 1.1σh. This initial homeostatic simulation ensured that subsequent G&R would be induced almost entirely at the location(s) where mass loss is introduced–that is, the modeled site(s) of injury.

Fig. 2
Deviation of stress relative to homeostatic stress (σkσh)h before (left) and after (right) 14 time steps of “accelerated” G&R simulation for the 45°(a)(b) and 135°(c)(d) collagen ...

While convergence of wall tension σ to homeostatic value can be obtained via adjustment of geometry or mass density distribution, convergence of wall shear stress τw can only be obtained through adjustment of vascular geometry. Because the initial homeostatic simulation should not introduce large geometric change, i.e., the original geometry should be maintained, wall shear stress regulation was not included in the initial homeostatic simulation.

3.2 Simulation for aneurysm expansion

An interesting application of the coupled G&R framework is to simulate aneurysm expansion in a region of vascular damage. Aneurysm expansion can be induced through constituent mass loss at specific locations, coupled with sufficiently low mass production rate of vascular constituents. Initial mass loss was introduced for the following three scenarios:

  1. Banded mass loss distribution
    Mk(x,0)={Mhk(x)(1-0.5fmr(cos(π(z-z0)2R0)+1)),z0-2R0zz0+2R0,Mhk(x),else.
  2. Point mass loss distribution
    Mk(x,0)={Mh(k)(x)(1-fmr),x-xcR0,Mhk(x)(1-fmr(2R0-x-xcR0)),R0<x-xc2R0,Mhk(x),else.
  3. Multiple point mass loss distribution
    Mk(x,0)={Mhk(x)(1-fmr),x-xc10.5R0,Mhk(x)(1-fmr(R0-x-xc10.5R0)),0.5R0<x-xc1R0,Mhk(x)(1-fmr),x-xc20.5R0,Mhk(x)(1-fmr(R0-x-xc20.5R0)),0.5R0<x-xc2R0,Mhk(x),else.
    Mhk(x) is the mass density generated by the initial homeostatic simulation. fmr is the maximum mass loss percentage set herein to 50%. R0 is the nominal value of the radius of the abdominal aorta, which was equal to 0.75 cm.

Figures 4(a)(b), 4(c)(d), and 4(e)(f) display the different initial mass loss distributions, along with the resulting vascular geometries after 110 G&R steps assuming Kσ = 0. (Below, we consider nontrivial values for Kσ.) Each G&R step was set to 0.1 of the collagen lifespan, so that 110 steps is equivalent to ≈800 days, assuming collagen has a life span around 72 days. It can be observed that significant expansion is only induced within or near the region of initial mass loss while the homeostatic state is effectively maintained in all other regions.

Fig. 4
Initial mass loss distribution(ratio between current and homeostatic mass density) and resulting aneurysm shapes(color denotes values of det(F)) after 110 G&R steps. Correspondence relations: (a)→(b), (c)→(d), (e)→(f)

COMSOL was used to solve the Navier Stokes equations. For illustrative purpose, a steady inflow was used with a resting infrarenal abdominal aortic flow rate equal to 15.2 cc/s. This approach considered time-averaged flow conditions since we are interested in the affect of the flow over time scales much longer than the cardiac cycle. Alternatively, we could apply a pulsatile inflow boundary condition and then perform time-averaging of the resulting wall shear stress and pressure fields. The outlets of the 3D fluid domain were at the level of the iliac arteries. To set the resistance values for each outlet, an equivalent resistance was first computed so that the mean pressure was 100 mmHg. This resistance was then distributed in parallel to the iliac arteries, resulting in R = 1.8269 ×109 Pa·s/m3, and R = 1.8285 × 109 Pa·s/m3, for the left and right iliac, respectively.

To explore the importance of feedback gains on vascular growth in our coupled framework, different values were considered for Kτ and Kσ. These results are summarized in Fig. 5, where the point mass loss shown in Fig. 4(c) was used to trigger G&R in all cases. All panels in this figure represent the respective field and geometry at the end of 110 G&R steps. Namely, columns are ordered by increasing expansion resulting from feedback gain choice (not evolution in time). The first column considers relatively large feedback gains: Kτ = 0.005 and Kσ = 0.005. In this case, the gains are large enough that even with the introduction of significant mass loss, expansion can be strongly suppressed through appropriate increase in wall thickness (mass density). After the 110 G&R steps, stress in both collagen constituents converged to the homeostatic value σh = 133kPa (see Fig. 6), and collagen turnover did not introduce any significant expansion.

Fig. 5
Results from coupled simulation with blood flow. The mass loss shown in Fig. 4(c) was used to trigger G&R in all cases. All panels represent the respective field and geometry at the end of 110 G&R steps, with each column representing different ...
Fig. 6
Convergence of constituent stress Kτ = 0.005, Kσ = 0.005 at the location of largest mass loss x = xc

We next considered four WSS feedback gains Kτ = 0.005, 0, −0.005, −0.01, whose results are respectively shown in columns 2–5 of Fig. 5. In each of these cases, we set Kσ = 0 to isolate the influence of WSS regulation on aneurysm progression, aortic pressure and blood flow. (Note that the case Kσ ≠ 0 and Kτ = 0 was considered previously in Fig. 4.) It can be observed that Kτ > 0 (corresponding to negative feedback) led to the lowest rate of expansion (Fig. 5, column 2). Positive feedback (corresponding to Kτ < 0) produced the most rapid expansion (Fig. 5, column 4 & 5). In all cases, it can be observed that WSS changed by roughly one order of magnitude within the aneurysm region, which contributed to G&R, whereas the change in transmural pressure varied very little (less than 0.1%).

It can be observed from Fig. 5 that both positive and negative values for Kτ led to expansion. Note that in the absence of any feedback (Kσ, Kτ = 0, Fig. 5, column 3), expansion occurs due to the normal turnover of collagen. Namely, the initial insult leads to immediate stretching of the existing collagen so that wall tension balances transmural pressure. This collagen is in an elevated stress state. Natural turnover of collagen produces new collagen whose prestress is lower than this elevated stress state, which causes these new fibers to stretch. This process drives a “basal expansion” due to the insult and natural collagen turnover. Kτ > 0 helps to suppress this basal expansion (column 2), whereas Kτ < 0 contributes to further expansion (columns 4, 5).

4 Discussion

We have presented a framework for coupling vascular growth and remodeling with blood flow simulation in a 3D patient-specific geometry using a constrained mixture theory for the vessel wall. Prior computational models [6, 18], coupling blood flow with G&R simulation have been applied to idealized geometry, or to a cylindrical segment of realistic geometry. Vascular G&R intrinsically occurs in all regions, and the purpose of this paper was to develop coupled simulation of G&R with blood flow in a more complete and anatomically realistic vascular model than had previously been considered.

One difficulty in applying G&R simulations to complex vascular domains derived from medical image data is defining appropriate initial conditions. When starting from a healthy model, it is reasonable to assume an initial homeostatic configuration. We demonstrated, as similarly shown in [6, 32] that an accelerated G&R stage can be performed to produce an approximate homeostatic configuration without significant geometrical alteration. As shown in Fig. 2, the relative stress deviation was nontrivial (>10%) in some regions after the initial accelerated G&R stage. The rate of convergence of the stress to the homeostatic value is proportional to the deviation. Therefore relatively long simulations times are needed to obtain very small deviations, whereas herein a relatively modest number of G&R steps were used. Indeed, our later simulation for Kσ = 0.005, Kτ = 0.005 verified that stresses do eventually converge to very near the homeostatic values after 110 G&R steps. Moreover, our simulations (Fig. 5) also show that the initial nontrivial deviations in stress from the homeostatic values do not result in significant unbalanced G&R for the time scales considered herein, and that significant expansion is almost entirely induced near or within the region where mass loss is introduced. Regardless, it might be more natural in future implementations to apply a homeostatic range instead of specifying a single homeostatic value.

Another challenge of applying G&R to complex vascular domains is consistently defining local axial and circumferential directions, with respect to which local anisotropic material properties are defined. In [31] this was addressed by defining local axial and circumferential directions using a 2-D parameterization of the vessel wall surface. However, such parameterization cannot be defined for domains that include multiple/branching vessels. To overcome this issue herein, local principal curvature directions on the surfaces were used to represent local axial and circumferential directions. Initial collagen directions (45° and 135°) were defined with respect to axial and circumferential directions in the reference configuration κ0 and later evolve with the mixture deformation.

In several prior coupled G&R modeling studies, aneurysm progression is triggered by prescribing initial insult as in [6]. Time dependent insult models have also been used [26] whereby AAA formation is initiated by prescribing a spatial and temporal degradation function for elastin. In more recent papers [27, 1, 17, 18], elastin degradation has been coupled to hemodynamics (mainly WSS) in order to initiate aneurysm progression. Herein, we used a simple insult model to trigger aneurysm progress since our focus was on application of coupled simulation to patient-specific geometry. However, more sophisticated initiation processes can be incorporated in this framework when applied to specific clinical applications.

In the coupled simulations with blood flow, the influence of G&R feedback gains Kσ and Kτ on the resulting aneurysm geometry and hemodynamics were studied. For large values of Kσ = 0.005 and Kτ = 0.005 (column 1, Fig. 5), growth and remodeling can compensate for relatively large initial mass loss so that the resulting geometry changed very little. Fig. 6 shows the convergence of the stress σk(t) to the homeostatic value σh at the center of mass loss, xc. We note that xc is the location of largest mass loss and generally corresponds to largest deformation; hence recovery of σk at xc generally indicates the convergence of stress in entire geometry.

We set Kσ = 0 to promote aneurysm expansion and then varied the value of Kτ to study the influence of the coupling with WSS regulation on the aneurysm progression. The role of WSS in aneurysm progression is complex and not well understood, therefore we considered cases for Kτ being positive, negative and zero. Indeed, prior publications have considered Kτ positive [23] and negative [6]. Mass loss will necessarily lead to some expansion when Kσ = 0. However, Kτ > 0 will compensate and reduce expansion rate due to the induced negative feedback. And conversely, expansion was exaggerated when Kτ < 0 due to the induced positive feedback.

The WSS plots in Fig. 5 show that the magnitude of WSS is very low in the aneurysm region (less than 0.05Pa), consistent with reported values in AAA [16]. We note that even for Kτ > 0 WSS does not converge to the homeostatic value. This may be because WSS regulation was turned off during the initial accelerated G&R stage or due to the fact that WSS is more sensitive to vascular geometry and hence more difficult to converge. Lastly, we note that it was verified that setting Kτ = 0 led to identical expansion to a case that did not consider coupling with fluid dynamics. From the pressure plots in Fig. 5, it can be observed that pressure change within the 3-D computational domain is small (less than 0.5%), and pressure changes even less at each location between the cases considered. Therefore, it may be reasonable to apply a uniform pressure for all cases instead of explicit coupling the pressure information from the fluid simulation.

The choice of material parameters describing the passive vascular response given in Table 1 were obtained from [32]. These parameters are based on a phenomenological model, which may included smooth muscle. Applying these parameters in the model herein, which neglects the passive response due to smooth muscle, may introduce additional uncertainty to these parameters. However, because smooth muscle is at least an order of magnitude less stiff than collagen, its contribution to the passive stress-strain behavior is expected to be small. Hence, we do not anticipate that such additional uncertainty in parameters will be significant, especially in comparison to the inherent uncertainty in these parameters.

5 Conclusions

Herein a computational framework to couple vascular growth and remodeling with blood flow simulation in a 3D patient-specific geometry was presented. We demonstrated that stress mediated regulation of wall tension and wall shear stress led to expected long-term response in vessel wall progression following mass loss for different feedback parameters. These results extend prior computational work on coupling G&R with hemodynamics simulation to more complex, subject-specific setting. Additionally, by coupling the time scales of hemodynamics and G&R simulation we are able to better understand the connections between local, short term hemodynamic factors and long term vascular change. Indeed, most image-based blood flow modeling has been concerned with understanding hemodynamic phenomena occurring on the time scale of the cardiac cycle for diagnosing current flow conditions. However, there is compelling need to develop simulation capabilities to predict blood flow conditions and vascular adaption over time scales of months to years, which is difficult, or impossible, to consider experimentally. This framework helps to bridge this gap.

Acknowledgments

We would like to thank Dr. Seungik Baek for helpful discussions during the development of this work. This work was supported in part by National Heart, Lung, and Blood Institute Grant 5R21-HL-108272.

Appendix

Algorithm for coupled simulation of G&R and hemodynamics

  1. Generate initial homeostatic state through accelerated G&R.
    1. Introduce initial mass loss to the homeostatic state.
    2. Simulate blood flow in the initial geometry and obtain initial WSS field τw(x, t).
    3. Set initial values for mass production rate mimk(t).
  2. Formulate the weak form for displacement u = (u, v, w) from the virtual work principle,
    δI=SδwdA-sPn·δxda=0,
    and solve using a nonlinear FEM solver.
    1. Generate stress measure field σk(x, t) based on deformation and constitutive relations
    2. Update mass production rate:
      mk(t)=M(t)M(0)(Kσ(σk(t)-σh)-Kτ(τw(t)-τwh)+fh)
  3. If ||mk(t) − mi||/||mi|| < tolerance (in this work, tolerance= 0.001),
    1. Set mimk(t)
    2. Set tt + Δt
    3. Update values for remaining fractions Qk(·) and qk(·)
    4. Go to Step (6)
    Else,
    1. Set mimk(t)
    2. Go to Step (3)
  4. If the geometric change is significant for fluid domain,
    1. Solve blood flow in the new geometry
    2. Update WSS field τw(x, t) and pressure P(x, t)
    3. Go to Step (3)
    Else, go to Step (3)

Footnotes

Conflicts of Interest

The authors do not have conflicts of interest relevant to this manuscript.

References

1. Aparício P, Mandaltsi A, Boamah J, Chen H, Selimovic A, Bratby M, Uberoi R, Ventikos Y, Watton PN. Modelling the influence of endothelial heterogeneity on the progression of arterial disease: Application to abdominal aortic aneurysm evolution. International Journal for Numerical Methods in Biomedical Engineering. 2014;30(5):563–586. [PubMed]
2. Baek S, Rajagopal KR, Humphrey JD. Competition between radial expansion and thickening in the enlargement of an intracranial saccular aneurysm. Journal of Elasticity. 2005;80(1–3):13–31.
3. Baek S, Rajagopal KR, Humphrey JD. A theoretical model of enlarging intracranial fusiform aneurysms. Journal of Biomechanical Engineering. 2006;128(1):142–149. [PubMed]
4. Burton AC. Relation of structure to function of the tissues of the wall of blood vessels. Physiol Rev. 1954;34(6) [PubMed]
5. Esmaily Moghadam M, Vignon-Clementel IE, Figliola R, Marsden AL. A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations. Journal of Computational Physics. 2013;244:63–79.
6. Figueroa AC, Baek S, Taylor CA, Humphrey JD. A computational framework for fluid–solid-growth modeling in cardiovascular simulations. Computer Methods in Applied Mechanics and Engineering. 2009;198(45):3583–3602. [PMC free article] [PubMed]
7. Finlay Helen M, Whittaker Peter, Canham Peter B. Collagen organization in the branching region of human brain arteries. Stroke. 1998;29(8):1595–1601. [PubMed]
8. Hariton I, Christian Gasser T, Holzapfel GA, et al. Stress-modulated collagen fiber remodeling in a human carotid bifurcation. Journal of theoretical Biology. 2007;248(3):460–470. [PubMed]
9. Holzapfel GA, Gasser TC, Stadler M. A structural model for the viscoelastic behavior of arterial walls: continuum formulation and finite element analysis. European Journal of Mechanics-A/Solids. 2002;21(3):441–463.
10. Humphrey JD, Holzapfel GA. Mechanics, mechanobiology, and modeling of human abdominal aorta and aneurysms. Journal of Biomechanics. 2012;45(5):805–814. [PMC free article] [PubMed]
11. Humphrey JD, Rajagopal KR. A constrained mixture model for growth and remodeling of soft tissues. Mathematical Models and Methods in Applied Sciences. 2002;12(03):407–430.
12. Malek A, Izumo S. Physiological fluid shear stress causes downregulation of endothelin-1 mrna in bovine aortic endothelium. American Journal of Physiology-Cell Physiology. 1992;263(2):C389–C396. [PubMed]
13. Niedermüller H, Skalicky M, Hofecker G, Kment A. Investigations on the kinetics of collagen-metabolism in young and old rats. Experimental Gerontology. 1977;12(5):159–168. [PubMed]
14. Rizvi MAD, Katwa L, Spadone DP, Myers PR. The effects of endothelin-1 on collagen type i and type iii synthesis in cultured porcine coronary artery vascular smooth muscle cells. Journal of Molecular and Cellular Cardiology. 1996;28(2):243–252. [PubMed]
15. Rizvi MAD, Myers PR. Nitric oxide modulates basal and endothelin-induced coronary artery vascular smooth muscle cell proliferation and collagen levels. Journal of Molecular and Cellular Cardiology. 1997;29(7):1779–1789. [PubMed]
16. Salsac A-V, Sparks SR, Chomaz J-M, Lasheras JC. Evolution of the wall shear stresses during the progressive enlargement of symmetric abdominal aortic aneurysms. Journal of Fluid Mechanics. 2006;560:19–51.
17. Selimovic A, Ventikos Y, Watton PN. Modelling the evolution of cerebral aneurysms: Biomechanics, mechanobiology and multiscale modelling. Procedia IU-TAM. 2014;10:396–409.
18. Sheidaei A, Hunley SC, Zeinali-Davarani S, Raguin LG, Baek S. Simulation of abdominal aortic aneurysm growth with updating hemodynamic loads using a realistic geometry. Medical Engineering & Physics. 2011;33(1):80–88. [PubMed]
19. Taber LA. A model for aortic growth based on fluid shear and fiber stresses. Journal of Biomechanical Engineering. 1998;120(3):348–354. [PubMed]
20. Uematsu M, Ohara Y, Navas JP, Nishida K, Murphy TJ, Alexander RW, Nerem RM, Harrison DG. Regulation of endothelial cell nitric oxide synthase mrna expression by shear stress. American Journal of Physiology-Cell Physiology. 1995;38(6):C1371. [PubMed]
21. Valentín A, Cardamone L, Baek S, Humphrey JD. Complementary vasoactivity and matrix remodelling in arterial adaptations to altered flow and pressure. Journal of The Royal Society Interface. 2009;6(32):293–306. [PMC free article] [PubMed]
22. Valentín A, Holzapfel GA. Constrained mixture models as tools for testing competing hypotheses in arterial biomechanics: A brief survey. Mechanics Research Communications. 2012;42:126–133. [PMC free article] [PubMed]
23. Valentín A, Humphrey JD. Evaluation of fundamental hypotheses underlying constrained mixture models of arterial growth and remodelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2009;367(1902):3585–3606. [PMC free article] [PubMed]
24. Valentín A, Humphrey JD, Holzapfel GA. A finite element-based constrained mixture implementation for arterial growth, remodeling, and adaptation: Theory and numerical verification. International Journal for Numerical Methods in Biomedical Engineering. 2013;29(8):822–849. [PMC free article] [PubMed]
25. Vande Geest JP, Sacks MS, Vorp DA. Age dependency of the biaxial biomechanical behavior of human abdominal aorta. Journal of Biomechanical Engineering. 2004;126(6):815–822. [PubMed]
26. Watton PN, Hill NA, Heil M. A mathematical model for the growth of the abdominal aortic aneurysm. Biomechanics and Modeling in Mechanobiology. 2004;3(2):98–113. [PubMed]
27. Watton PN, Raberger NB, Holzapfel GA, Ventikos Y. Coupling the hemodynamic environment to the evolution of cerebral aneurysms: computational framework and numerical examples. Journal of Biomechanical Engineering. 2009;131(10):101003. [PubMed]
28. Wilson JS, Baek S, Humphrey JD. Importance of initial aortic properties on the evolving regional anisotropy, stiffness and wall thickness of human abdominal aortic aneurysms. Journal of The Royal Society Interface. 2012;9(74):2047–2058. [PMC free article] [PubMed]
29. Wilson JS, Baek S, Humphrey JD. Parametric study of effects of collagen turnover on the natural history of abdominal aortic aneurysms. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science. 2013;469(2150):20120556. [PMC free article] [PubMed]
30. Wilson NM, Ortiz AK, Johnson AB. The vascular model repository: A public resource of medical imaging data and blood flow simulation results. Journal of Medical Devices. 2013;7(4):040923. [PMC free article] [PubMed]
31. Zeinali-Davarani S, Baek S. Medical image-based simulation of abdominal aortic aneurysm growth. Mechanics Research Communications. 2012;42:107–117.
32. Zeinali-Davarani S, Sheidaei A, Baek S. A finite element model of stress-mediated vascular adaptation: application to abdominal aortic aneurysms. Computer Methods in Biomechanics and Biomedical Engineering. 2011;14(9):803–817. [PubMed]