Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC4497867

Formats

Article sections

Authors

Related links

Ann Biomed Eng. Author manuscript; available in PMC 2016 July 1.

Published in final edited form as:

Published online 2015 March 3. doi: 10.1007/s10439-015-1287-6

PMCID: PMC4497867

NIHMSID: NIHMS669011

Mechanical Engineering, University of California, Berkeley, CA, United States

Correspondence: Shawn C. Shadden, 5126 Etcheverry Hall, University of California, Berkeley, California 94720-1740, Email: ude.yelekreb@neddahs

The publisher's final edited version of this article is available at Ann Biomed Eng

See other articles in PMC that cite the published article.

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.

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.

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} *κ*(*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 **G*** ^{k}*(

$${\mathbf{F}}_{n(\tau )}^{k}(t)=\mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){\mathbf{G}}^{k}(\tau ),$$

(1)

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

$${\mathbf{C}}_{n(\tau )}^{k}(t)={\mathbf{F}}_{n(\tau )}^{k}{(t)}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathbf{F}}_{n(\tau )}^{k}(t).$$

(2)

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

The unit vector in the direction of collagen family *k* at time *t* is denoted **e*** ^{k}*(

$${\mathbf{e}}_{n(\tau )}^{k}=\frac{{\mathbf{G}}^{k}{(\tau )}^{-1}{\mathbf{e}}^{k}(\tau )}{\Vert {\mathbf{G}}^{k}{(\tau )}^{-1}{\mathbf{e}}^{k}(\tau )\Vert}.$$

(3)

With **e*** ^{k}*(

$${\mathbf{G}}^{k}(\tau )={G}^{k}\phantom{\rule{0.16667em}{0ex}}{\mathbf{e}}^{k}(\tau )\otimes {\mathbf{e}}_{n(\tau )}^{k},$$

(4)

where *G ^{k}* is the stretch ratio of newly produced collagen fiber from the natural configuration
${\kappa}_{n(\tau )}^{k}$ to the mixture configuration

$$\begin{array}{l}{\lambda}_{n(\tau )}^{k}(t)=\sqrt{{\mathbf{e}}_{n(\tau )}^{k}\xb7{\mathbf{C}}_{n(\tau )}^{k}(t){\mathbf{e}}_{n(\tau )}^{k}}\\ =\Vert \mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){\mathbf{G}}^{k}(\tau ){\mathbf{e}}_{n(\tau )}^{k}\Vert \\ =\Vert \mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){G}_{h}^{c}\left({\mathbf{e}}^{k}(\tau )\otimes {\mathbf{e}}_{n(\tau )}^{k}\right)\phantom{\rule{0.16667em}{0ex}}{\mathbf{e}}_{n(\tau )}^{k}\Vert \\ ={G}_{h}^{c}\Vert \mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert \\ ={G}_{h}^{c}\frac{\Vert \mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert /\Vert {\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert}{\Vert {\mathbf{e}}^{k}(\tau )\Vert /\Vert {\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert}\\ ={G}_{h}^{c}\frac{{\lambda}^{k}(t)}{{\lambda}^{k}(\tau )},\end{array}$$

(5)

where *λ ^{k}*(

$${\lambda}^{k}(t)=\frac{\Vert \mathbf{F}(t){\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert}{\Vert {\mathbf{F}}^{-1}(\tau ){\mathbf{e}}^{k}(\tau )\Vert}.$$

(6)

For elastin, we assume the mapping from the natural configuration to the mixture reference configuration is **G ^{e}**, and therefore the mapping from elastin’s natural configuration to the current mixture configuration is given by

$${\mathbf{F}}_{n}^{e}(t)=\mathbf{F}(t){\mathbf{G}}^{e}.$$

(7)

In this model we take ${\mathbf{G}}^{\mathbf{e}}=\text{diag}[{G}_{1}^{e},{G}_{2}^{e},{\scriptstyle \frac{1}{{G}_{1}^{e}\phantom{\rule{0.16667em}{0ex}}{G}_{2}^{e}}}]$, 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.

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 *M ^{k}*(

$${M}^{k}(t)={M}^{k}(0){Q}^{k}(t)+\underset{0}{\overset{t}{\int}}{m}^{k}(\tau ){q}^{k}(t-\tau )d\tau .$$

(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. *Q ^{k}*(

For elastin,

$${M}^{e}(t)={M}^{e}(0){Q}^{e}(t),$$

(9)

where *Q ^{e}*(

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

$${q}^{k}(t)=\{\begin{array}{ll}1,\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}0\le t<{t}_{1}\hfill \\ {\scriptstyle \frac{1}{2}}\left\{cos\phantom{\rule{0.16667em}{0ex}}\left({\scriptstyle \frac{\pi}{{t}_{2}-{t}_{1}}}(t-{t}_{1})\right)+1\right\},\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{t}_{1}\le t\le {t}_{2}\hfill \\ 0,\hfill & \text{if}\phantom{\rule{0.16667em}{0ex}}{t}_{2}<t,\hfill \end{array}$$

(10)

where *t*_{2} is the lifespan of the constituent. For collagen we use *t*_{1}*/t*_{2} = 0.2, and *t*_{2} = 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

$${Q}^{k}(t)=1-\frac{{\int}_{0}^{t}{q}^{k}(\tau )d\tau}{{\int}_{0}^{{t}_{2}}{q}^{k}(\tau )d\tau}=1-\frac{2{\int}_{0}^{t}{q}^{k}(\tau )d\tau}{{t}_{1}+{t}_{2}}.$$

(11)

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

$${W}^{k}\left({I}_{n(\tau )}^{k}\right)=\frac{{c}_{2}}{4{c}_{3}}\left\{exp\phantom{\rule{0.16667em}{0ex}}\left[{c}_{3}\phantom{\rule{0.16667em}{0ex}}{\left({I}_{n(\tau )}^{k}-1\right)}^{2}\right]-1\right\},$$

(12)

where

$${I}_{n(\tau )}^{k}={\left({\lambda}_{n(\tau )}^{k}(t)\right)}^{2}.$$

(13)

Let

$${\mathbf{C}}_{n(\tau )}^{k}(t)={\mathbf{F}}_{n(\tau )}^{k}{(t)}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathbf{F}}_{n(\tau )}^{k}(t)$$

(14)

denote the right Cauchy-Green deformation tensor of the *k*th family of collagen produced at time *τ*, whose current direction is denoted by **e*** ^{k}*(

$${W}^{e}({\mathbf{C}}_{n}^{e}(t))=\frac{{c}_{1}}{2}({I}_{1}-3),$$

(15)

where

$${I}_{1}(t)=\text{tr}({\mathbf{C}}_{n}^{e}(t)).$$

(16)

${\mathbf{C}}_{n}^{e}(t)={\mathbf{F}}_{n}^{e}{(t)}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathbf{F}}_{n}^{e}(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

$$\begin{array}{l}w(t)={w}^{e}(t)+\sum _{k}{w}^{k}(t)\\ ={M}^{e}(0){Q}^{e}(t){W}^{e}({\mathbf{C}}_{n}^{e}(t))+\sum _{k}\left\{{M}^{k}(0){Q}^{k}(t){W}^{k}\left({I}_{n(0)}^{k}\right)+\underset{0}{\overset{t}{\int}}{m}^{k}(\tau ){q}^{k}(t-\tau ){W}^{k}\left({I}_{n(\tau )}^{k}\right)\right\}\end{array}$$

(17)

where *w ^{e}*(

$$\delta I={\int}_{S}\delta w\phantom{\rule{0.16667em}{0ex}}dA-{\int}_{s}P\phantom{\rule{0.16667em}{0ex}}\mathbf{n}\xb7\delta \mathbf{x}\phantom{\rule{0.16667em}{0ex}}da=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.

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
${\mathbf{e}}_{0}^{k}$ 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 **e*** ^{k}*(

$${\mathbf{e}}^{k}(t)=\mathbf{F}(t){\mathbf{e}}_{0}^{k}.$$

(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.

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

$${\rho}_{b}\phantom{\rule{0.16667em}{0ex}}\left(\frac{\partial \mathbf{v}(\mathbf{x},t)}{\partial t}+\mathbf{v}(\mathbf{x},\mathbf{t})\xb7\nabla \mathbf{v}(\mathbf{x},\mathbf{t})\right)=-\nabla p(\mathbf{x},t)+\mu {\nabla}^{2}\mathbf{v}(\mathbf{x},\mathbf{t}),$$

(20)

$$\nabla \xb7\mathbf{v}(\mathbf{x},\mathbf{t})=0.$$

(21)

The blood density *ρ _{b}* and viscosity

$${P}_{o}(t)=R{Q}_{o}(t)$$

(22)

The flow rate *Q _{o}* at the respective outlets was obtained from the 3D domain, and using Eq. (22) the pressure

It is important to note that the time scales for G&R (weeks) and blood flow simulation (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.

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

$$\begin{array}{l}\mathbf{T}(t)=\frac{1}{J(t)}\mathbf{F}(t)\frac{\partial w}{\partial {I}_{1}(t)}\frac{\partial {I}_{1}(t)}{\partial \mathbf{F}(t)}+\frac{1}{J(t)}\mathbf{F}(t)\sum _{k}\frac{\partial w}{\partial {I}_{n(\tau )}^{k}(t)}\frac{\partial {I}_{n(\tau )}^{k}(t)}{\partial \mathbf{F}(t)}+{\mathbf{T}}_{\text{active}}\\ =\frac{2}{J(t)}\frac{\partial {w}^{e}}{\partial {I}_{1}(t)}\mathbf{B}(t)+\frac{2}{J(t)}\sum _{k}\frac{\partial {w}^{k}}{\partial {I}_{n(\tau )}^{k}}\frac{\partial {I}_{n(\tau )}^{k}}{\partial {I}^{k}(t)}{\mathbf{e}}^{k}\otimes {\mathbf{e}}^{k}+{\mathbf{T}}_{\text{active}},\end{array}$$

(23)

where **B** = **FF**^{} 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 **T**_{active} 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

$${\sigma}^{k}=\frac{{\mathbf{e}}^{k}\xb7{\mathbf{Te}}^{k}}{h},$$

(24)

which is the stress in the **e*** ^{k}* direction. The thickness of the vessel wall was calculated as

$$h(t)=\frac{M(t)}{J\rho},$$

(25)

where *M* (*t*) = Σ* _{k} M^{k}*(

In G&R theory, the vascular homeostatic state is recovered through stress-mediated feedback. The mass production rate of the *k*th family of collagen is assumed to depend linearly on the deviation of wall tension *σ ^{k}* with respect to the homeostatic value

$${m}^{k}(t)=\frac{M(t)}{M(0)}\left({K}_{\sigma}\phantom{\rule{0.16667em}{0ex}}\left({\sigma}^{k}(t)-{\sigma}^{h}\right)-{K}_{\tau}\phantom{\rule{0.16667em}{0ex}}\left({\tau}_{w}(t)-{\tau}_{w}^{h}\right)+{f}_{h}^{\stackrel{\sim}{k}}\right),$$

(26)

where
${\stackrel{\sim}{f}}_{h}^{k}$ is the basal value of mass production rate for collagen family *k*,

$${\stackrel{\sim}{f}}_{h}^{k}=\frac{{M}^{k}(0)}{{\int}_{0}^{\infty}{q}^{k}(\tau )d\tau},$$

(27)

which balances the degradation rate of collagen in the homeostatic state. Note that *τ _{w}*(

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].

Local collagen directions calculated based on principal curvature directions at the aortic bifurcation: (a) 45° collagen family, (b) 135° collagen family.

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

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.

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:

- Banded mass loss distribution$${M}^{k}(\mathbf{x},0)=\{\begin{array}{ll}{M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}\left(1-0.5{f}_{mr}\phantom{\rule{0.16667em}{0ex}}\left(cos\phantom{\rule{0.16667em}{0ex}}\left(\frac{\pi (z-{z}_{0})}{2{R}_{0}}\right)+1\right)\right),\hfill & {z}_{0}-2{R}_{0}\le z\le {z}_{0}+2{R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x}),\hfill & \text{else}.\hfill \end{array}$$
- Point mass loss distribution$${M}^{k}(\mathbf{x},0)=\{\begin{array}{ll}{M}_{h}^{(k)}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}(1-{f}_{mr}),\hfill & \mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}}\mid \phantom{\rule{0.16667em}{0ex}}\le {R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}\left(1-{f}_{mr}\phantom{\rule{0.16667em}{0ex}}\left(\frac{2{R}_{0}-\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}}\mid}{{R}_{0}}\right)\right),\hfill & {R}_{0}<\phantom{\rule{0.16667em}{0ex}}\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}}\mid \phantom{\rule{0.16667em}{0ex}}\le 2{R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x}),\hfill & \text{else}.\hfill \end{array}$$
- Multiple point mass loss distribution$${M}^{k}(\mathbf{x},0)=\{\begin{array}{ll}{M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}(1-{f}_{mr}),\hfill & \mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{1}}\mid \phantom{\rule{0.16667em}{0ex}}\le 0.5{R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}\left(1-{f}_{mr}\phantom{\rule{0.16667em}{0ex}}\left(\frac{{R}_{0}-\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{1}}\mid}{0.5{R}_{0}}\right)\right),\hfill & 0.5{R}_{0}<\phantom{\rule{0.16667em}{0ex}}\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{1}}\mid \phantom{\rule{0.16667em}{0ex}}\le {R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}(1-{f}_{mr}),\hfill & \mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{2}}\mid \phantom{\rule{0.16667em}{0ex}}\le 0.5{R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x})\phantom{\rule{0.16667em}{0ex}}\left(1-{f}_{mr}\phantom{\rule{0.16667em}{0ex}}\left(\frac{{R}_{0}-\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{2}}\mid}{0.5{R}_{0}}\right)\right),\hfill & 0.5{R}_{0}<\phantom{\rule{0.16667em}{0ex}}\mid \mathbf{x}-{\mathbf{x}}_{\mathbf{c}\mathbf{2}}\mid \phantom{\rule{0.16667em}{0ex}}\le {R}_{0},\hfill \\ {M}_{h}^{k}(\mathbf{x}),\hfill & \text{else}.\hfill \end{array}$$${M}_{h}^{k}(\mathbf{x})$ is the mass density generated by the initial homeostatic simulation.
*f*is the maximum mass loss percentage set herein to 50%._{mr}*R*_{0}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

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 *×*10^{9} Pa·s/m^{3}, and *R* = 1.8285 *×* 10^{9} Pa·s/m^{3}, 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

Convergence of constituent stress *K*_{τ} = 0.005, *K*_{σ} = 0.005 at the location of largest mass loss **x** = **x**_{c}

We next considered four WSS feedback gains *K _{τ}* = 0.005

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 (

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

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

We set *K _{σ}* = 0 to promote aneurysm expansion and then varied the value of

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

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.

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.

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.

Algorithm for coupled simulation of G&R and hemodynamics

- Generate initial homeostatic state through accelerated G&R.
- Introduce initial mass loss to the homeostatic state.
- Simulate blood flow in the initial geometry and obtain initial WSS field
*τ*(_{w}**x***, t*). - Set initial values for mass production rate
*m*→_{i}*m*(^{k}*t*).

- Formulate the weak form for displacement
**u**= (*u, v, w*) from the virtual work principle,$$\delta I={\int}_{S}\delta \mathit{wdA}-{\int}_{s}P\mathbf{n}\xb7\delta \mathbf{x}da=0,$$and solve using a nonlinear FEM solver. - Generate stress measure field
*σ*(^{k}**x***, t*) based on deformation and constitutive relations - Update mass production rate:$${m}^{k}(t)=\frac{M(t)}{M(0)}\left({K}_{\sigma}\phantom{\rule{0.16667em}{0ex}}\left({\sigma}^{k}(t)-{\sigma}^{h}\right)-{K}_{\tau}\phantom{\rule{0.16667em}{0ex}}\left({\tau}_{w}(t)-{\tau}_{w}^{h}\right)+{\stackrel{\sim}{f}}_{h}\right)$$

- If ||
*m*(^{k}*t*) −*m*||_{i}*/*||*m*||_{i}*<*tolerance (in this work, tolerance= 0.001),- Set
*m*←_{i}*m*(^{k}*t*) - Set
*t*←*t*+ Δ*t* - Update values for remaining fractions
*Q*(·) and^{k}*q*(·)^{k} - Go to Step (6)

Else,- Set
*m*←_{i}*m*(^{k}*t*) - Go to Step (3)

- If the geometric change is significant for fluid domain,
- Solve blood flow in the new geometry
- Update WSS field
*τ*(_{w}**x***, t*) and pressure*P*(**x***, t*) - Go to Step (3)

Else, go to Step (3)

**Conflicts of Interest**

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

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |