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

**|**HHS Author Manuscripts**|**PMC4778980

Formats

Article sections

- Abstract
- 1 Introduction
- 2 Methods
- 3 Results
- 4 Discussion
- 5 Limitations and Future Work
- 6 Conclusions
- References

Authors

Related links

Theor Comput Fluid Dyn. Author manuscript; available in PMC 2017 April 1.

Published in final edited form as:

Published online 2015 December 19. doi: 10.1007/s00162-015-0374-5

PMCID: PMC4778980

NIHMSID: NIHMS746555

Vittoria Flamini, Department of Mechanical and Aerospace Engineering, New York University Tandon School of Engineering, Brooklyn, New York, USA;

The publisher's final edited version of this article is available at Theor Comput Fluid Dyn

See other articles in PMC that cite the published article.

It has long been recognized that aortic root elasticity helps to ensure efficient aortic valve closure, but our understanding of the functional importance of the elasticity and geometry of the aortic root continues to e-volve as increasingly detailed in vivo imaging data become available. Herein, we describe a fluid-structure interaction model of the aortic root, including the aortic valve leaflets, the sinsuses of Valsalva, the aortic annulus, and the sinotubular junction, that employs a version of Peskin’s immersed boundary (IB) method with a finite element (FE) description of the structural elasticity. As in earlier work, we use a fiber-based model of the valve leaflets, but this study extends earlier IB models of the aortic root by employing an incompressible hyperelastic model of the mechanics of the sinuses and ascending aorta using a constitutive law fit to experimental data from human aortic root tissue. In vivo pressure loading is accounted for by a backward displacement method that determines the unloaded configurations of the root model. Our model yields realistic cardiac output at physiological pressures, with low transvalvular pressure differences during forward flow, minimal regurgitation during valve closure, and realistic pressure loads when the valve is closed during diastole. Further, results from high-resolution computations indicate that although the detailed leaflet and root kinematics show some grid sensitivity, our IB model of the aortic root nonetheless produces essentially grid-converged flow rates and pressures at practical grid spacings for the high-Reynolds number flows of the aortic root. These results thereby clarify minimum grid resolutions required by such models when used as stand-alone models of the aortic valve as well as when used to provide models of the outflow valves in models of left ventricular fluid dynamics.

The aortic root consists of the three semilunar aortic valve cusps, the three sinuses of Valsalva, which are bulbous cavities positioned behind each valve leaflet, the aortic annulus, where the aortic valve leaflets attach to the aorta, and the sinotubular junction, where the sinuses merge into the ascending aorta. The healthy aortic valve acts to ensure unidirectional flow of oxygenated blood from the heart to the tissues of the body, including to the heart itself via the coronary circulation. Diseases of the aortic valve can result in stenosis or regurgitation, and severe aortic valve disease is treated by repairing or by replacing the diseased valve with either a mechanical or a bioprosthetic valve [11]. Approximately 50,000 such procedures are performed each year in the United States [6, 26, 67]. Continuing advances in noninvasive in vivo imaging of blood flow and tissue deformation, mechanical tests specifically targeted to biological tissues, and computer modeling and simulation enable integrative studies of the dynamics and mechanics of the aortic root [41]. Such work has the potential to improve both medical devices and surgical procedures used to treat patients with valvular heart disease, and also clinical approaches to patient risk assessment and treatment planning.

Aortic compliance, and specifically the compliance of the aortic sinuses, has a fundamental role in the function of the aortic root. The sinuses act as reservoirs, storing blood during systole and then releasing it in diastole to facilitate flow in the coronary arteries [71]. The compliance of the sinuses also helps to ensure the proper closure of the aortic valve [4]. Further, the shape of the aortic root causes blood to circulate in the sinuses, generating vortices that act both on the valve leaflets and on the sinus walls. These vortices generate forces that facilitate effective valve closure and act as a regulatory mechanism that synchronizes closure [4]. Although this mechanism of efficient valve closure was first postulated by Leonardo da Vinci, there remains some debate in the surgical community about whether it is important to recreate the anatomic geometry of the native aortic root in aortic valve and root replacement procedures [12].

More recently, our understanding of the functional importance of the elasticity and anatomical geometry of the aortic root has evolved as increasingly detailed data have been acquired on the in vivo deformations that occur during the cardiac cycle. For instance, it has been shown that the aortic root undergoes a multi-modal series of conformational changes even before the leaflets open [18]. It has been argued that these deformations act to reduce shear stresses on the valve leaflets and thereby prolong the life of the native valve leaflets [12]. Annular and commissural flexibility may be a key component in this interaction, and annular flexibility is lost with the implantation of a stented artificial valve, possibly reducing the lifetime of bioprosthetic valve leaflets. The annular commissures also participate in load sharing, which reduces peak stresses on the aortic cusps immediately after valve closure. Devices such as the Medtronic 3f Aortic Bioprosthesis attempt to account for both the annular and the commissural flexibility of the native root, but despite the theoretical advantages offered by such designs, more traditional devices with rigid annuluses still dominate clinical practice.

A challenge in modeling and simulating the mechanical response of arteries, including the aorta, is that the arterial walls are continuously subjected to substantial intraluminal pressures in their in vivo state. Thus, arterial geometries that are derived from in vivo imaging generally correspond to a loaded configuration. In addition, as first demonstrated by Fung [27] and by Vaishnav and Vossoughi [72], even the unloaded configuration of the artery is not stress free. Instead, because of tissue growth and remodeling, the arterial wall includes residual stresses and is subject to axial tethering [10, 28, 41]. These factors all conspire to make determining the unloaded arterial geometry challenging. Indeed, both the deformation due to intraluminal pressure loading and the residual stresses cannot be measured directly in vivo, and therefore these quantities must be estimated. Residual stresses may be determined from ex vivo experiments, whereas the zero-pressure configuration of blood vessels can be determined by solving the inverse elastostatic problem.

Different approaches have been developed to solve the inverse elastostatic problem for complex arterial geometries derived from medical imaging data. Most of these strategies can be grouped in two broad categories. One group of approaches focuses on retrieving the initial deformation field and the initial stress field of the artery while keeping the imaging-derived geometry unaltered. An example of this type of approach is the modified updated Lagrangian formulation (MULF) introduced by Gee et al. [31], which is an incremental prestressing method that recovers the equilibrium configuration rather than the zero-pressure geometry. A second group of approaches determines an unloaded configuration of the vessel that, when subjected to intraluminal pressure loading, will inflate to match the imaging-derived geometry. Among the first to develop a stationary method to solve the inverse elastostatic problem were Govindjee and Mihalic [32]. This method was extended by Lu et al. [46] to arteries, but this approach does not appear to have been widely used in practice, possibly because it requires a modification in the finite element (FE) solution scheme that renders this method somewhat difficult to implement within commercial FE analysis software. In comparison, the backward incremental method developed by de Putter et al. [20] requires updating the deformation gradient and thereby can be implemented with fewer modifications to the finite element code. Perhaps the most straightforward approach to retrieving the zero-pressure configuration is the backward displacement method proposed by Bols et al. [5], which iteratively modifies the coordinates of the reference configuration rather than the deformation gradient tensor. We use the backward displacement method in this work to recover the unloaded configuration of an idealized model of the aortic root.

While the methods to obtain an estimation of the zero-pressure configuration of a blood vessel can be applied to any geometry, determining residual stresses is a local, vessel-specific task that presently requires harvesting and then destroying the blood vessel. Computational approaches to accounting for residual stresses are usually based on an assumed ‘open’ configuration, which is subsequently ‘closed,’ thereby reversing the process that releases residual stresses [2, 15]. This approach is difficult to apply to noncylindrical geometries and, in particular, cannot be readily applied to imaging-derived geometries [15]. A more generic approach to including residual stress in arteries is described in Alastrué et al. [2] and is based on a Kröner-Lee decomposition of the deformation gradient tensor. For simplicity, we do not consider residual stresses in the present work, although we do account for initial strains resulting from in vivo pressure loading.

Most the methods surveyed above have been developed and applied to study aneurysmal arteries [31, 46], although applications to normal vessels as well as other pathophysiological conditions are possible. Because arterial wall stress distributions are considered good predictors of aneurysm rupture, improvements in the accuracy of computed wall stress distributions are expected to facilitate improvements in aneurysm rupture prediction. In contrast, although the important role of aortic root compliance in valve closure is well documented [16, 17], many dynamic analyses of the aortic valve tend to focus on the mechanical behavior of the valve leaflets and often model the sinuses and the ascending aorta as rigid [34, 38, 74, 78] or as linear elastic materials [14, 49], but experimental data [3] show that the sinuses and the ascending aorta exhibit a nonlinear hyperelastic behavior. Examples of nonlinear hyperelastic constitutive models used in fluid-structure interaction simulations of the aortic root can be found in studies by De Hart et al. [19] and by Weinberg and Mofrad [75], which both adopt an arbitrary Lagrangian-Eulerian (ALE) approach. The application of ALE fluid-structure interaction schemes to modeling aortic valve dynamics is somewhat limited by the remeshing required by such methods as the structure deforms [68]. Attempts to tackle this limitation have included restricting the problem to a two-dimensional analysis [23] and the implementation of complex remeshing algorithms in three-dimensional analyses [13].

An alternative approach to modeling interactions between the blood and the aortic valve and root is offered by Peskin’s immersed boundary (IB) method [58], which was introduced to simulate cardiac valve dynamics [57] and was subsequently extended to model fluid-structure interaction in the heart and its valves [34–36,38,44,47,48,51–53,60]. This IB approach to fluid-structure interaction is to describe structural stresses and deformations in Lagrangian form, and to describe the momentum, viscosity, and incompressibility of the coupled fluid-structure system in Eulerian form. Integral transforms with Dirac delta function kernels mediate interaction between Lagrangian and Eulerian variables. When discretized, the IB formulation of the equations of motion replaces these singular kernel functions with regularized kernels that are designed to ensure conservation of physical quantities such as force and torque when converting between Lagrangian and Eulerian forms. This numerical approach allows the Lagrangian structural mesh to overlay the background Eulerian grid in an arbitrary manner, thereby avoiding the need to deploy dynamic body-fitted grids. In addition, because the immersed structures move according to a common interpolated velocity field, the IB method offers an implicit contact model that prevents the valve leaflets from interpenetrating, even when subjected to substantial diastolic pressure loads [34, 38]. Related sharp-interface IB methods have also been developed and used by Borazjani, Ge, and Sotiropoulos to simulate valvular dynamics for models of both rigid and flexible aortic valve prostheses [7, 8]. Previous work demonstrated that three-dimensional IB models of the aortic valve can yield physiological cardiac output at realistic pressures [34, 38]. However, we are aware of no previous study that demonstrates that Peskin’s IB method yields reasonably well-resolved simulation results for flexible aortic valve models in three spatial dimensions at practical grid spacings.

The primary aim of this study is to develop a new fluid-structure interaction model of the aortic root that substantially extends earlier IB models of aortic valve dynamics [34,38] by including descriptions of the elasticity, distensibility, and dynamics of the aortic sinuses and the ascending aorta. Herein, the aortic sinuses and proximal ascending aorta are modeled as an incompressible, isotropic, hyperelastic material with an exponential neo-Hookean strain-energy functional [21, 41] that we fit to experimental data obtained by Azadani et al. [3] from human aortic root tissue samples. A separate finite element analysis is employed to estimate the unloaded aortic geometry using a method based on the backward displacement method of Bols et al. [5]. Our fluid-structure interaction simulations employ hybrid models of the aortic root that use a fiber-based description of the thin aortic valve leaflets, as done in previous studies [34, 38], along with a finite element-based description [39] of the comparatively thick walls of the aortic sinuses and proximal ascending aorta.

In our dynamic analysis, we pace the aortic root model to an essentially periodic steady state using a prescribed, periodic left-ventricular pressure waveform. A human Windkessel model fit to clinical data [69] provides downstream loading for the aortic root, and a pressure waveform, taken from the clinical data set used to parameterize this reduced circulation model [55], serves to drive flow through the model. In the overall model, flow rates are not prescribed, but rather emerge from the computation. Physiological cardiac output is obtained at physiological driving and loading pressures, with low transvalvular pressure differences during forward flow, minimal regurgitation during valve closure, and realistic transvalvular pressure loads when the valve is closed during diastole.

The aortic root model considered in this work captures the complex interactions between the flow, the thin flexible aortic valve cusps, and the deformable walls of the aortic sinuses and the ascending aorta. It employs an idealized anatomical geometry, in which the initial and reference configurations of the aortic sinuses and valve cusps are assumed to exhibit three-fold symmetry, but can be extended to more realistic anatomic geometries. This model is fully three dimensional in the sense that there are no symmetry conditions imposed on the subsequent fluid dynamics or structural kinematics. Using this somewhat idealized model, we perform a refinement study to demonstrate that the present methods are able to yield essentially grid-converged flow rates and pressures at practical grid spacings. Specifically, under grid refinement, only relatively small differences are observed in stroke volume, maximum flow velocity, and vessel distensibility and stress distribution. Relatively large differences remain in the fine details of the flow in the vicinity of the valve cusps, and in the details of the kinematics of the valve cusps and the aortic wall. Higher spatial resolution is likely needed to resolve fully the dynamics of the valve during systole. Nonetheless, a key contribution of this study is that it demonstrates that bulk flow properties are reasonably resolved by IB models of aortic valve dynamics at practical spatial resolutions. These results are practically useful for determining minimum grid resolutions when using IB models of the outflow valves.

The immersed boundary formulation used herein describes fluid-structure
systems in which an elastic structure is immersed in a viscous incompressible
fluid. It employs a Lagrangian description of the structural deformations and
the stresses generated by those deformations, and an Eulerian description of the
momentum, incompressibility, and viscosity of the coupled fluid-structure
system. In the present model, we employ a fiber-based description of the thin
valve leaflets, and we describe the aortic wall as an incompressible
hyperelastic solid. Let **q** *U*
^{2} indicate Lagrangian curvilinear coordinates attached
to the valve leaflets, with **q** = (*q, r*), let
**X** *V* ^{3}
indicate the Lagrangian material coordinate system of the aortic wall, with
**X** = (*X, Y, Z*), and let **x**
Ω indicate the Eulerian physical coordinates of the physical
domain, with **x** = (*x, y, z*). The physical position
of fiber point **q** at time *t* is given by
**ϕ**(**q**, *t*)
Ω, and the physical position of aortic wall material point
**X** at time *t* is given by
**χ**(**X**, *t*)
Ω. We use Ω^{l}(*t*) Ω
to indicate the (codimension 1) physical region occupied by the valve leaflets
at time *t*, Ω^{w}(*t*)
Ω to indicate the (codimension 0) physical region occupied by the
solid-body model of the vessel wall, and
Ω^{f}(*t*) = Ω \
Ω^{w}(*t*) to indicate the physical region
occupied by the fluid at time *t*.

The equations of motion for the coupled fluid-structure system are:

$$\rho \frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t}(\mathbf{x},t)=-\nabla p(\mathbf{x},t)+\mu {\nabla}^{2}\mathbf{u}(\mathbf{x},t)+\mathbf{f}(\mathbf{x},t)+\mathbf{g}(\mathbf{x},t),$$

(1)

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

(2)

$$\mathbf{f}(\mathbf{x},t)={\int}_{U}\mathbf{F}(\mathbf{q},t)\delta (\mathbf{x}-\mathit{\varphi}(\mathbf{q},t))\mathrm{d}\mathbf{q},$$

(3)

$$\mathbf{g}(\mathbf{x},t)={\int}_{V}{\nabla}_{\mathbf{X}}\xb7\mathbb{P}(\mathbf{X},t)\delta (\mathbf{x}-\mathit{\chi}(\mathbf{X},t))\mathrm{d}\mathbf{X}\phantom{\rule{0ex}{0ex}}-{\int}_{\partial V}\mathbb{P}(\mathbf{X},t)\mathbf{N}(\mathbf{X})\delta (\mathbf{x}-\mathit{\chi}(\mathbf{X},t))\mathrm{d}A(\mathbf{X}),$$

(4)

$$\frac{\partial \mathit{\varphi}}{\partial t}(\mathbf{q},t)={\int}_{\Omega}\mathbf{u}(\mathbf{x},t)\delta (\mathbf{x}-\mathit{\varphi}(\mathbf{q},t))\mathrm{d}\mathbf{x},$$

(5)

$$\frac{\partial \mathit{\chi}}{\partial t}(\mathbf{X},t)={\int}_{\Omega}\mathbf{u}(\mathbf{x},t)\delta (\mathbf{x}-\mathit{\chi}(\mathbf{X},t))\mathrm{d}\mathbf{x},$$

(6)

in which ρ is the mass density, μ is the dynamic
viscosity, **u**(**x**, *t*) is the Eulerian
velocity field of the fluid-structure system, *p*(**x**,
*t*) is the Eulerian pressure field of the fluid-structure
system, **f** (**x**, *t*) is the Eulerian
elastic force density generated by deformations to the fiber model of the valve
leaflets, **g**(**x**, *t*) is the Eulerian
elastic force density generated by deformations to the solid-body model of the
vessel wall, **F**(**q**, *t*) is the
Lagrangian elastic force density of the fiber model, (**X**,
*t*) is the first Piola-Kirchhoff elastic stress tensor of
the solid model, **N**(**X**) is the outward unit normal along
*V*, δ(**x**) =
δ(*x*) δ(*y*)
δ(*z*) is the three-dimensional Dirac delta function,
and $\frac{\mathrm{D}(\xb7)}{\mathrm{D}t}=\frac{\partial (\xb7)}{\partial t}+\mathbf{u}\xb7\nabla (\xb7)$ is the material derivative. In the present work, we assume
uniform mass density ρ and dynamic viscosity μ. These
assumptions are not essential, however, and versions of the IB method have been
developed that permit the use of spatially varying mass densities [25, 42, 43, 54, 62, 79] and viscosities [25, 62].

Eqs. (1) and (2) are the Eulerian incompressible
Navier-Stokes equations. Here, the momentum equation (1) is augmented by two Eulerian body force
densities. The first of these, **f** (**x**,
*t*), is determined by eq. (3) to be the Eulerian force density that is equivalent
to the Lagrangian fiber force density **F**(**q**,
*t*). The second body force in (1), **g**(**x**, *t*), is
determined by eq. (4) to be the
Eulerian force density equivalent to the Lagrangian description of the forces
generated by the solid-body model of the vessel wall, which are expressed in
terms of (**X**, *t*). Notice that
**f** (**x**, *t*) is a singular force
density supported on Ω^{l}(*t*). By contrast,
**g**(**x**, *t*) is supported on
Ω^{w}(*t*). Away from
Ω^{w}(*t*),
**g**(**x**, *t*) is nonsingular, although
**g**(*x, t*) is singular along
Ω^{w}(*t*) where
**N** ≠ **0**; see eq. (4).

Eqs. (5) and (6) determine the motions of the
immersed structures from the Eulerian material velocity field
**u**(**x**, *t*). Because of the presence
of viscosity throughout the fluid-solid system, **u**(**x**,
*t*) is continuous, and thus (5) and (6)
are equivalent to

$$\frac{\partial \mathit{\varphi}}{\partial t}(\mathbf{q},t)=\mathbf{u}(\mathit{\varphi}(\mathbf{q},t),t),$$

(7)

$$\frac{\partial \mathit{\chi}}{\partial t}(\mathbf{X},t)=\mathbf{u}(\mathit{\chi}(\mathbf{X},t),t).$$

(8)

Because · **u**(**x**,
*t*) = 0, the immersed solid is automatically treated as
incompressible in this formulation. Specifically, at least in the continuum
formulation, there is no need to include terms in the elastic strain-energy
functional to penalize compressible deformations. Notice that this formulation
describes the inertia of both the fluid and the solid, and the immersed solid
bodies are not assumed or required to be in (quasi-static) equilibrium.

In practice, we use a standard *C*^{0} finite
element method to approximate the deformations and stresses of the vessel wall.
To do so, it is useful to adopt a weak formulation for the structural equations,
so that

$$\mathbf{g}(\mathbf{x},t)={\int}_{V}\mathbf{G}(\mathbf{X},t)\delta (\mathbf{x}-\mathit{\chi}(\mathbf{X},t))\mathrm{d}\mathbf{X},$$

(9)

$${\int}_{V}\mathbf{G}(\mathbf{X},t)\xb7\mathbf{V}(\mathbf{X})\mathrm{d}\mathbf{X}=-{\int}_{V}\mathbb{P}(\mathbf{X},t):{\nabla}_{\mathbf{X}}\mathbf{V}(\mathbf{X})\mathrm{d}\mathbf{X},\forall \mathbf{V}(\mathbf{X}),$$

(10)

$$\frac{\partial \mathit{\chi}}{\partial t}(\mathbf{X},t)=\mathbf{U}(\mathbf{X},t),$$

(11)

$${\int}_{V}\mathbf{U}(\mathbf{X},t)\xb7\mathbf{V}(\mathbf{X})\mathrm{d}\mathbf{X}={\int}_{V}\mathbf{u}(\mathbf{\chi}(\mathbf{X},t),t)\xb7\mathbf{V}(\mathbf{X})\mathrm{d}\mathbf{X},$$

(12)

in which **G**(**X**, *t*) is
the Lagrangian elastic force density of the aortic wall,
**U**(**X**, *t*) is the Lagrangian
velocity field of the aortic wall, and **V**(**X**) is an
arbitrary Lagrangian test function that is not assumed to vanish on
*V*. In the continuum setting, eqs. (4) and (9)–(10) are equivalent definitions
for **g**(**x**, *t*); however, these
formulations lead to different numerical schemes when discretized, and only
eqs. (9)–(10) lead directly to a standard
nodal finite element structural discretization. Further, in the continuum
equations, eqs.
(11)–(12) are
equivalent to eq. (6) and to
eq. (8). Further details of
the numerical methods used in this work are described in a brief Appendix.

As in earlier work [34, 38], the leaflet geometry is determined by the mathematical theory of the fiber architecture of aortic valve leaflets developed by Peskin and McQueen [59]. This theoretical model of valve geometry and fiber architecture derives the shape of the leaflets from their function, which is to support a uniform pressure load during diastole, when the valve is closed. Each valve leaflet is composed of two families of fibers that are orthogonal to each other. One family of fibers runs from commissure to commissure, and the other runs from the bottom scalloped edge of the aortic sinus to the free edge of the valve leaflet. The leaflets are defined in a closed and loaded configuration; see fig. 1. In this configuration, the radius of each leaflet, measured from the tip of the valve to the end of the belly region, is 1.45 cm, and the coapting portion of each leaflet is 0.97 cm tall. The scale of the leaflets is based on measurements of human aortic roots [61,70]. The height of the coapting portion of the leaflet is chosen to be similar to that of the real valve while enabling the model valve to support a physiological pressure load when closed.

The mechanical behavior of the leaflets is described by a
strain-energy functional *E* =
*E*[**ϕ**(·,
*t*)] described previously [34, 38]. Briefly, let the
curvilinear coordinates **q** = (*q, r*) be chosen
so that a fixed value of *q* labels an individual fiber, and
so that *r* runs along the fibers. Consequently,
(**ϕ**/*r*)/‖**ϕ**/*r*‖
is the unit fiber tangent vector. The total elastic energy
*E* is the sum of a stretching energy
*E*_{s} and a bending energy
*E*_{b},

$$E={E}_{\mathrm{s}}+{E}_{\mathrm{b}},$$

(13)

$${E}_{\mathrm{s}}={\int}_{U}{\mathcal{E}}_{\mathrm{s}}\left(\right|\frac{\partial \mathit{\varphi}}{\partial r}|;\mathbf{q})\mathrm{d}\mathbf{q},$$

(14)

$${E}_{\mathrm{b}}=\frac{1}{2}{\int}_{U}{c}_{\mathrm{b}}(\mathbf{q}){|\frac{{\partial}^{2}\mathit{\varphi}}{\partial {r}^{2}}-\frac{{\partial}^{2}\overline{\mathit{\varphi}}}{\partial {r}^{2}}|}^{2}\mathrm{d}\mathbf{q}.$$

(15)

Eq. (14) accounts for the
total stretching energy associated with the fibers, in which
_{s} is a spatially inhomogeneous local stretching
energy with a quadratic length-tension relationship, and eq. (15) accounts for the
total bending energy associated with the fibers, in which
*c*_{b} is a spatially inhomogeneous bending
stiffness and $\overline{\mathit{\varphi}}$ is the reference configuration, which is taken to be the
initial configuration. Because the valve leaflets are modeled as thin
elastic surfaces, the bending resistant energy allows the model to account
for the thickness of the real valve leaflets. Larger values of
*c*_{b} model a thick, stenotic valve, whereas
smaller values of *c*_{b} model a thin, flexible
valve. The resulting Lagrangian body force density is the Fréchet
derivative of *E*.

Valve leaflet parameters are empirically determined for a fixed leaflet discretization by choosing the leaflet stiffness so that the leaflet supports a diastolic load of approximately 80 mmHg with 10% strain in the family of fibers running from commissure to commissure. This leads to a commissural fiber stiffness of 7.5e6 dyne/cm. We further assume that these commissural fibers are a factor of 10 stiffer than the radial fibers that run orthogonal to the commissural fibers [64]. The bending stiffness is chosen to be approximately the smallest value that yields valve leaflets that successfully coapt at the end of systole. These material parameters are chosen only to yield a functional valve, and are not represented as corresponding to those of an actual valve cusp. Nonetheless, this model of the valve leaflets does account for key structures of actual valve cusps, including the prominent collagen fibers that run from commissure to commissure and the 10-to-1 anisotropy of real valve leaflets [64]. See refs. 38 and 34 for further details, and see also the discussion of the limitations of this model in sec. 5. In preliminary studies, a sensitivity analysis was conducted for the bending stiffness of the aortic leaflets, and it was found that variations on the order of ± 25% in the bending stiffness have little effect on the physiological predictions of the models (not shown).

An idealized model of the aortic root is constructed, as in previous work [34,38]. The dimensions of this model are based on measurements by Swanson and Clark [70] of human aortic roots collected at autopsy inflated to a diastolic pressure of 80 mmHg, and the geometry of the model is based on measurements by Reul et al. [61] from angiograms. The diameter of the aortic portion of the model is 3 cm, whereas the diameter of the sinus region, measured from a commissura to the center of a sinus, is 3.5 cm. The overall length of the model is of 10 cm, and the distance between the annulus and the aortic flow outlet is 7.75 cm. The thickness is constant throughout the model and is 2 mm [66]. In our simulations, we employ a semi-rigid model of the left-ventricular outflow tract while allowing for a fully flexible description of the aortic sinuses and the ascending aorta. The aortic annulus, which is the scalloped line of attachment between the valve leaflets and the aortic sinuses, was considered to be the lower boundary of the sinuses; see fig. 1. The aortic model used in this study differs from previous work [34,38] in that here we treat the aortic sinuses and ascending aorta as an incompressible hyperelastic material, whereas in earlier studies it was treated as a rigid structure.

Because the geometry of the valve leaflets is defined in a closed
and loaded configuration, it is also necessary to specify the geometry of
the vessel wall in a loaded configuration, and to compute the corresponding
unloaded configuration. This configuration will depend on both the
prescribed loaded configuration and the constitutive model.^{1} Because we use a constitutive
model fit to data obtained from aortic tissue samples with a
“J-shaped” stress-strain relationship (see fig. 2), accounting for initial
deformations is also needed to obtain realistic deformations at systemic
pressure loads. Specifically, if we did not account for these initial
deformations, the sinuses and ascending aorta would experience
unrealistically large deformations under even diastolic pressures. To
determine the unloaded configuration, we employ an iterative method
described by Bols et al. [5] detailed
in sec. 2.3.3.

We describe the elasticity of the aortic sinuses and ascending aorta
using a hyperelastic constitutive model fit to biaxial tensile test data
collected by Azadani et al. [3] from
tissue samples from human aortic sinuses. Because the data reported by
Azadani et al. indicate a nearly isotropic material response (see fig. 2), in this work we model the
material response of the aortic sinuses and ascending aorta using an
isotropic strain-energy functional *W* with an exponential
stress-strain relation,

$$W=\frac{c}{2b}[\text{exp}(b({I}_{1}-3))-1],$$

(16)

in which *I*_{1} =
*I*_{1}() = tr() is the first
invariant of the right Cauchy-Green strain tensor =
^{T} , and =
**χ**/**X** is the
deformation gradient tensor associated with the deformation mapping
**χ** : (*V, t*) Ω.
Material parameters *c* and *b* were obtained
by least-squares fits to experimental data using MATLAB (Mathworks, Inc.,
Natick, MA, USA). The first Piola-Kirchhoff stress tensor is
determined from *W* via

$$\mathbb{P}=\frac{\partial W}{\partial \mathbb{F}}-{p}_{\mathrm{s}}{\mathbb{F}}^{-T}=c\text{exp}(b({I}_{1}-3))\mathbb{F}-{p}_{\mathrm{s}}{\mathbb{F}}^{-T},$$

(17)

in which *p*_{s} is a structural
pressure-like term that is chosen to improve the accuracy of the stress
predictions of the method [29]. As in
earlier work [29], we compute
*p*_{s} via

$${p}_{\mathrm{s}}=c\text{exp}(b({I}_{1}-3))-{\beta}_{\mathrm{s}}\text{log}({I}_{3}),$$

(18)

in which *I*_{3} =
*I*_{3}() = det() is the third
invariant of and β_{s} = 2.5e4 kPa. Thus,
= **0** for = , and the structural
model provides an energetic penalty for any compressible deformations.

To determine the unloaded configuration of the vessel wall, we use an iterative backward displacement method described by Bols et al. [5] implemented using a custom MATLAB script that interfaces with the ABAQUS (Dassault Systèmes Simulia Corp., Providence, RI, USA) finite element analysis software. Because the constitutive model (16) is not provided by ABAQUS, this model was implemented using a UHYPER subroutine.

Let **x**_{m} denote positions of
the finite element mesh nodes in the *prescribed* loaded
configuration, let ${\mathbf{X}}_{m}^{i}$ denote the corresponding nodal positions in the
*computed* reference configuration after
*i* iterations of the backward displacement algorithm,
and let ${\mathit{\chi}}_{m}^{i}={\mathit{\chi}}_{m}[{\mathbf{X}}^{i}]$ denote nodal positions in the *computed*
loaded configuration when **X**^{i} is
used as the unloaded reference configuration. Notice that
**x**_{m}, ${\mathbf{X}}_{m}^{i}$, and ${\mathit{\chi}}_{m}^{i}$ are all defined at the nodes of the finite element
mesh.

The algorithm is straightforward: First, we initialize
**X**^{0} := **x** (i.e., we use the
deformed coordinates as an initial guess for the unloaded configuration).
Then, given **X**^{i},
**X**^{i+1} is determined by

$${\mathbf{X}}^{i+1}\u2254{\mathbf{X}}^{i}+(\mathbf{x}-{\mathit{\chi}}^{i}).$$

(19)

Notice that if we define the displacement from the computed
reference configuration to the deformed coordinates at step
*i* by **D**^{i} =
**χ**^{i} −
**X**^{i}, then (19) is equivalent to

$${\mathbf{X}}^{i+1}\u2254\mathbf{x}-{\mathbf{D}}^{i}.$$

(20)

Thus, the reference coordinates are determined via a
*backward* displacement from the deformed configuration.
If this process converges,
**χ**^{i} →
**x** and we obtain reference coordinates **X** such
that **χ**[**X**] = **x**.

The convergence of this iterative process is assessed in terms of
the discrete *L*^{2} norm of the residual
**x** −
**χ**^{i}, i.e.

$${r}^{i}={\Vert \mathbf{x}-\mathit{\chi}[{\mathbf{X}}^{i}]\Vert}_{2}={\left(\sum _{m}{\Vert {\mathbf{x}}_{m}-{\mathit{\chi}}_{m}^{i}\Vert}_{2}^{2}\right)}^{1/2}.$$

(21)

The routine was judged to have converged once the residual was less than 0.1% of the average vessel radius and the change in the residual between two consecutive iterations was less than 0.015% of the average radius. In general, obtaining convergence may require the use of damping, or of gradually increasing the loading pressure; however, in practice, we have found this algorithm to be robust, and neither damping nor incremental loading were required for the analyses performed herein.

A left-ventricular pressure waveform adapted from human clinical data of Murgo et al. [55] is prescribed at the upstream inlet of the left-ventricular outflow tract to drive flow through the model aortic root, and downstream loading conditions are provided by a three-element Windkessel model by Stergiopolus et al. [69] fit to the clinical data of Murgo et al. [55]. As in earlier work [34,38], zero pressure boundary conditions are imposed on the boundary of the fluid-filled region exterior to the aortic root model. Coupling between the detailed IB model and the reduced circulation models is also done as in earlier work [34,38]. Along the proximal and distal ends of the vessel, zero-displacement boundary conditions are imposed.

Because the normal traction boundary conditions are used at both the
inlet and the outlet are subject to the development of flow instabilities [24], a small amount of flow stabilization
is applied in the fluid region within a distance of approximately
2*h* from the boundary. Specifically, the net flow rate
through the upstream or downstream boundary is first determined, and then any
local flow that is opposite from the direction of net flow is penalized by
including additional forcing that acts to avoid localized regions of flow
reversal on the boundary. The same basic approach is used both for inflow and
outflow, and it permits the boundary conditions to switch from inflow to outflow
conditions based on the flow within the interior of the detailed model, and on
the pressure differences across that model. This stabilization scheme is similar
in its physical mechanism to that described by Esmaily Moghadam [24], but here we employ interior forcing
rather than directly modifying the formulation of the physical boundary
condition.

In our dynamic fluid-structure interaction simulations, the
computational domain is taken to be a 10 cm × 10 cm × 10 cm
cubic region discretized using a three-level locally refined grid. For the
purposes of a mesh convergence study, two different fine grid spacings are used:
a coarser one corresponding to *h* = 0.78125 mm, which is
equivalent to a uniform 128 × 128 × 128 Eulerian discretization,
and a finer one corresponding to *h* = 0.390625 mm, which is
equivalent to a uniform 256 × 256 × 256 Eulerian discretization.
The curvilinear mesh used to discretize the valve leaflets has a physical grid
spacing of approximately *h*/4 for the coarser Eulerian
discretization, and to *h*/2 for the finer discretization. The
volumetric mesh used to discretize the vessel wall uses 23,400 8-node
(trilinear) hexahedral elements for the aortic root model, which results in an
average mesh aspect ratio of 1.93 ± 0.25. A preliminary convergence
study in ABAQUS was performed to ensure mesh convergence of the structural model
for both the backward displacement method and the forward FE analysis. We use
the standard four-point delta function of Peskin [58] for the coarser simulations, and a broadened 8-point
version of this kernel function for the higher resolution simulations, so that
the physical extent of the regularized delta functions is the same for both
spatial resolutions.^{2}
Consequently, the valve leaflets have the same effective thickness with respect
to the fluid dynamics part of the simulation for both spatial resolutions. We
remark that the structural response of the leaflets is primarily determined by
the structural model and not by the effective leaflet thickness
“seen” by the fluid model.

A uniform time step size of Δ*t* = 9.94369e-6 s
was used with the coarser Eulerian discretization, whereas a uniform time step
size of Δ*t* = 1.40625e-5 was used for the finer Eulerian
discretization. These time step sizes are empirically determined to be within a
factor of $\sqrt{2}$ of the largest time step that satisfies the stability
restriction of our explicit time integration method. CFL numbers are on the
order of 0.03 during systole and 0.003 during diastole. Thus, the time step size
is generally substantially smaller than the maximum stable time step size
allowed by our fluid solver, for which the CFL number is order 1. For this
model, the stability is primarily determined by the explicit coupling between
the solid and the Eulerian momentum equation. Asymptotically, we expect that for
this model, Δ*t* ~ *h*^{4}
because of the presence of bending-resistant elastic elements in the model valve
leaflets.

The constitutive model (16) used to describe the mechanical response of the aortic sinuses
and ascending aorta is fit to biaxial tensile tests of human aortic root tissue
samples reported by Azadani et al. [3],
yielding model parameter values of *c* = 12.8 kPa and
*b* = 6.9; see fig. 2.
Recall that these material parameters are used in the aortic sinuses and
ascending aorta, whereas the outflow tract is treated as an essentially rigid
structure.

The unloaded configuration is determined by assuming that the loaded configuration corresponds to a pressure load of 80 mmHg. In fig. 3, panel (a) shows the prescribed loaded geometry, and panel (b) shows the computed zero-pressure configuration. Fig. 4 shows the convergence history of the backward displacement method. When the computed zero-pressure geometry is subjected to a 80 mmHg pressure load, the deformed geometry agrees with the prescribed geometry to within approximately 15 µm; see fig. 5. The discrepancy between the prescribed loaded geometry and the computed loaded geometry is greatest at the commissures and smallest along the straight portion of the vessel.

A dynamic fluid-structure interaction analysis is performed that
includes four cardiac cycles, including an initialization cycle and three cycles
in which the model has attained periodic steady state. (Data are reported only
for the beats in which the model has reached steady state.) Fig. 6 shows the prescribed left-ventricular driving
pressure along with the computed aortic loading pressure determined by the
coupled Windkessel model and the computed flow rate through the symmetric aortic
root. We emphasize that the flow rate is *not* imposed in the
model, but rather is predicted by the fluid-structure interaction analysis.

Results of a dynamic fluid-structure interaction analysis of the aortic
root over the final three simulated cardiac cycles using a coarser Cartesian
grid (effective grid resolution of *N* = 128) and a finer
Cartesian grid (effective grid resolution of **...**

For the simulations using the coarser Cartesian spacing,
*h* = 0.78125 mm, stroke volume is 96.2 ml in average, peak
flow rate is 591.5 ml/s, and cardiac output is 6.5 l/min. For the finer
Cartesian grid spacing, *h* = 0.3906 mm, stroke volume is 100.2
ml in average, peak flow rate is 643.8 ml/s, and cardiac output is 6.8 l/min.
Both sets of simulation results are within 15% of the hemodynamic
parameters of the subject-specific clinical data used to determine the driving
pressure and the Windkessel model parameters. Specifically, the driving pressure
and the Windkessel model parameters used in this work were based on clinical
measurements from a subject with stroke volume of 100 ml, peak flow rate of 560
ml/s, and cardiac output of 6.8 l/min, as reported by Murgo et al. [55]. The simulations also capture the
incisura in the flow profile corresponding to the backward flow generated by the
leaflets. The backward flow volume is 1.77 ± 0.11 ml, or 1.8% of
the stroke volume, for the coarser grid spacing, and 1.52 ± 0.059 ml, or
1.5% of the stroke volume, for the finer grid spacing. The net forward
flow volume (i.e., into the aorta) following leaflet coaptation is 0.02
± 0.01 ml. This forward flow results from interactions between the
elastic aortic root and the dynamic downstream loading pressure. As the
downstream pressure falls during diastole, the load on the aortic wall and
leaflets decreases and allows the aortic root to gradually push fluid into the
distal aorta.

The model captures the opening and closing dynamics of the valve at both Cartesian grid spacings. Fig. 7 shows the leaflets being pushed apart along with the formation of the vortices that help the leaflets to close at the end of systole. Fig. 8 shows that although there is some net flow towards the left ventricle as the valve closes, once the valve closes, no further flow is allowed between the aorta and the left ventricle.

Valve leaflet configurations during valve opening shown at equally
spaced times during the fourth simulated cardiac cycle for (a) the coarser
Cartesian grid (effective grid resolution of *N* = 128) and (b)
the finer Cartesian grid (effective grid resolution **...**

Fig. 9 shows the distribution of
the structural displacements in the aortic root, and fig. 10 shows displacement as a function of time at four
selected locations: the center of the sinus; the sinotubular junction; the
commissure; and the ascending aorta. Notice that after the first cardiac cycle,
the displacements are essentially periodic in time. Cartesian grid resolution
has a minimal effect on the range of displacements except for the displacements
measured at the commissure. The largest amplitude is observed in the ascending
aorta (0.58 ± 0.05 mm) and at the sinotubular junction
(0.52±0.02 mm). A measure of aortic stiffness is the radial pulsation,
i.e., the ratio of the difference between the systolic radius and the diastolic
radius to the average radius [56], which
in the present analysis was approximately 3.9% for the coarser Cartesian
grid spacing and 2.9% for the finer Cartesian grid spacing. Notice that
radial pulsation is a highly sensitive measure of deformation, and the different
values obtained for the coarse and fine simulations indicate a relatively small
difference in the displacements of approximately 0.175 mm, which is less than
0.5*h* on the highest resolution Cartesian grid used in this
study. These values are similar to the value of 2.5% reported in refs.
41 and 56 for the ascending aorta. The distribution of the
displacements was also in good agreement with literature data, with the
displacements in the aortic root greater in the sinotubular junction, and at the
commissures than in the sinuses, as reported also by Lansac [45].

Displacements of the deformable aortic root model at four locations: (a)
the straight portion of the aorta; (b) just above one of the commissures; (c)
the sinotubular junction; and (d) the middle of the sinus.

The distribution of the maximum principal stress in the aortic root follows a distribution similar to that of the displacement, with the maximum stresses concentrated in the sinotubular junction, just above the sinuses, as shown in fig. 11. Stress in the aortic root follows the pulsatile waveform of the pressure; see fig. 12. Both Cartesian grid spacings yield similar stress distributions; compare fig. 11 panels (a) and (b), but fig. 12 shows some differences along the sinotubular junction and within the sinuses. The maximum principal stress in the sinuses is 20% of that observed in the rest of the aortic root. The distribution, and the range of the maximum principal stress predicted by our model, is in agreement with other analyses available in literature, in particular with the 220 kPa peak of the maximum principal stress in the sinotubular junction reported by Conti et al. [14].

In this work, previous IB models of aortic valve dynamics [34,38] are substantially extended by incorporating a deformable model of the aortic sinuses and ascending aorta. This extended model is thereby able to more closely replicate the complex in vivo dynamics of the real aortic root. As such, it promises to serve as a powerful model for simulations of aortic root behavior and valve dynamics in mechanobiology studies or in the design, testing, or selection of medical devices. The dynamic analysis of the aortic valve included four complete cardiac cycles, including a single initialization cycle, which are shown to be sufficient to yield essentially periodic results. As in previous IB models of aortic valve dynamics [34, 38], key hemodynamic parameters, including aortic pressure, maximum flow rate, cardiac output, and stroke volume, are in reasonable agreement with physiological ranges [22, 76] and are also in good agreement with the patient-specific data [55] used to fit the Windkessel model [69] employed in this study. Large oscillations in the flow rate, as well as in the deformations, indicate the reverberations of the aortic valve leaflets upon closure, and the resonances of the elastic model of the sinuses and ascending aorta; i.e., these are the so-called heart sounds that are generated by the aortic valve. Specifically, these dynamics are not the result of the explicit time coupling scheme used in the numerical method. Further, grid convergence results demonstrate for the first time that the IB method is able to yield essentially grid-resolved bulk flow parameters and vessel deformations at practical grid spacings.

Many strategies have been developed to determine the zero-pressure geometry of blood vessels [1, 5, 20, 31, 65, 73]. Indeed, medical images of arteries always provide the vessel geometry in a loaded configuration. It is known, however, that using the imaged geometry as the reference geometry of the vessel underestimates the amount of strain present in the blood vessel [20]. Among the approaches developed to estimate the zero-pressure geometry of blood vessels, we chose an iterative backward displacement approach, as described by Bols et al. [5] and by Sellier [65]. Iterative approaches based on geometric considerations are straightforward to implement and can be applied to a variety of problems and solvers because they require minimal implementation effort. In this work, the backward displacement analysis was implemented in ABAQUS and converged to a zero-pressure geometry that, when loaded, yielded a deformed configuration that is in good agreement with the initially prescribed loaded geometry; see figs. 3–5. However, the deflated geometry is quite different from the configuration assumed by a deflated aortic root harvested from a body; see fig. 3. This is a drawback of these methods, as shown also in a recent study by Wittek et al. [77], in which the unloaded configuration of the aorta is different from an unloaded aorta harvested from a body. This mismatch could be the result of many factors, such as the choice of the boundary conditions or the need for adding the residual stress to balance the level of stress within the aortic wall [28]. We expect that the addition of residual stress estimates, planned in future work, will reduce the mismatch between the deflated geometry computed and the deflated aortic root observed in vivo.

Interestingly, the addition of a compliant aortic wall did not substantially alter the global hemodynamic parameters previously determined by a noncompliant IB aortic root model [34]; see fig. 13, which shows that the peak pressure attained by the noncompliant model is slightly higher than that obtained in the compliant model, whereas the flow rates are similar in both models. This finding might suggest that our vessel wall model is overly stiff despite our use of a constitutive model fit to biaxial tensile test data from human aortic tissues samples. Another reason for the reduced motility of the aortic root model is the fact that the outflow tract of the left ventricle is fixed. In physiological conditions, due to ventricular contraction and relaxation, the aortic root moves upwards and sideways with every heartbeat. In this work, we constrained this motion, thus restricting the dynamics of the aortic root model. This constraint reduced greatly the accelerations the aortic root is subjected to as compared to the physiological conditions.

In contrast, we observed that using a compliant aorta affected more clearly the dynamics of valve leaflet closure. Figs. 7(b) and 8(b) show the flow inversion and the distension in the sinus, which collects fluid during systolic ejection, and the elastic recoil of the sinus, which pushes down on the leaflet during closure. The addition of this mechanism made the model substantially more robust. Specifically, in previous work [34] we found that the rigid aortic root model required the valve leaflet properties to be very carefully tuned in order for the model to remain competent throughout the cardiac cycle. This is in clear contrast to the present compliant aortic root model, which yields a competent valve for a relatively broad range of leaflet parameters (data not shown). A compliant root provides elastic recoil to “push” the leaflets closed and also makes it easier for the nearly-closed leaflets to deform to coapt fully and to prevent regurgitation. This is a trait that makes our model more relevant physiologically, as aortic valves are continuously remodeled and keep functioning in vivo for varying mechanical properties and configurations. Future work will address the impact that geometrical characteristics of the ascending aorta, such as length, curvature, and tortuosity, have on the FSI simulations. We expect that these changes will not substantially alter the bulk flow parameters calculated in this study, but we anticipate that changes in the shape of the ascending aorta will affect blood velocity distributions and will generate flow pattern qualitatively closer to those seen in vivo.

Similar model predictions were obtained for both Cartesian grid spacings used in this study. We observed hemodynamic values that were within 5% of the clinical values reported by Murgo et al. [55] except for the peak flow rate, which was within 5% of the values reported by Murgo et al. on the coarser Cartesian grid and within 15% on the finer grid. Cardiac output and stroke volume were both in very good agreement with the clinical data, within 4% on the coarser Cartesian grid, and within 1% on the finer grid. In addition, the pulsatile radius, a quantity used to measure aortic stiffness, was also shown to be in relatively good agreement with literature values in both cases [41].

Work is underway to address some of the limitations of these models. As described previously, our present models do not account for the curved geometry of the aortic root, or for the residual stresses present in real arterial vessels. We specifically aim to incorporate a description of residual stresses and realistic medical imaging-derived anatomic geometries into our IB models. Further, the present models use an isotropic model of the aortic sinuses and ascending aorta. Although the experimental data that were used to fit this model show an essentially isotropic material response, it is known that real aortic tissues have a well-defined collagen fiber structure with an anisotropic material response. We also aim to incorporate a more realistic hyperelastic model of the aortic wall that accounts for families of collagen fibers [30] as well as experimentally constrained models of the elasticity of the aortic valve leaflets [50]. Further, it is now well known that the valve leaflets are in fact layered structures [63]. These details are not accounted for in the present model, and we anticipate that including them in future work will require an important extension of the methods presented in this study. Specifically, because the valve leaflets are thin elastic structures, we expect that using realistic hyperelastic continuum models to describe the mechanics of the valve leaflets will require adopting a nonlinear shell formulation. To our knowledge, continuum shell models have not yet been widely used within the framework of the IB method, and we anticipate additional computational research will be required to implement these types of formulations into our IB simulation framework.

Another limitation of this work is the long run times required by these simulations. Because of the small time step sizes required by the present algorithm, each cardiac cycle required 50,000 or more time steps, requiring on the order of one day per cycle on 32–64 cores of a modern high-performance computing (HPC) cluster for the coarse simulations and on the order of one week per cycle for the fine-grid simulations. Higher-resolution simulations using the present methods and implementation are computationally infeasible at this time. We are currently working to develop effective multigrid solvers for stable implicit time stepping schemes for the IB method [40]. These methods promise to allow for substantially larger time step sizes as well as shorter run times at current grid resolutions. They could also enable higher resolution simulations that promise to resolve the fine-scale details eliminate the remaining grid-sensitivity observed in the kinematics of the valve leaflets and the distensibility of the aortic root.

This work has presented a new fluid-structure interaction model of the aortic root that extends earlier IB models of the aortic valve by incorporating a finite-strain model of the aortic root that uses a constitutive model fit to tensile test data obtained from human aortic root tissue specimens. This model captures both the complex fluid dynamics of the flow and the finite deformations of the vessel wall, and it yields results that are in good agreement with the clinical data that were used to fit the circulation model used in this study, as well as to physiological data in the research literature. Finally, a grid convergence study demonstrates that nearly grid converged results are obtained at practical spatial resolutions, although yet higher spatial resolution, or higher accuracy at present resolutions, is needed to obtain fully resolved simulation results.

This work has been funded in part by American Heart Association award 10SDG4320049, National Institutes of Health awards GM071558 and HL117063, and National Science Foundation awards DMS 1016554/1460368 (to NYU School of Medicine and UNC-Chapel Hill) and ACI 1047734/1460334 (to NYU School of Medicine and UNC-Chapel Hill). VF was supported in part by the Tandon School of Engineering of New York University. Computations were performed at New York University using computer facilities funded in large part by a generous donation by St. Jude Medical, Inc.

We summarize the key discretization methods used in this study in this appendix. We use a locally-refined staggered-grid discretization of the Eulerian equations along with a finite difference-based discretization of the fiber model of the valve leaflets and a finite element-based description of the continuum vessel wall model. Our approach is similar to the spatially adaptive IB scheme described previously [34], except that here we also employ a finite element-based description of the aortic wall, using an approach introduced by Griffith and Luo [39].

The physical domain Ω is discretized using a
locally-refined Cartesian grid, but for simplicity we describe only the
uniform-grid version of this method; details on the adaptively refined
discretization are provided by ref. 34. Let (*i, j, k*) index the cells of the
Cartesian grid, and let ${\mathbf{x}}_{i,j,k}=((i+\frac{1}{2})h,(j+\frac{1}{2})h,(k+\frac{1}{2})h)$ indicate the position of the center of grid cell
(*i, j, k*), in which *h* is the
Cartesian grid meshwidth and Δ**x** =
*h*^{3} is the Cartesian grid cell volume.
The Eulerian velocity field **u** = (*u*,
υ, *w*) is approximated at the center of each
face of the Cartesian grid cells in terms of the velocity component that
is normal to that face, so that *u* is approximated at
the locations ${\mathbf{x}}_{i-\frac{1}{2},j,k}$, υ is approximated at the locations ${\mathbf{x}}_{i,j-\frac{1}{2},k}$, and *w* is approximated at ${\mathbf{x}}_{i,j,k-\frac{1}{2}}$. The Eulerian force densities **f** and
**g** are approximated in the same staggered-grid fashion.
The Eulerian pressure *p* is approximated at the centers
of the Cartesian grid cells.

The deformations and forces associated with the valve leaflets
are approximated on a fiber-aligned curvilinear mesh. Let
*l* index the nodes of this mesh, and let
**ϕ**_{l} and
**F**_{l} indicate the current
position and Lagrangian force density of node *l*,
respectively, and let
Δ**q**_{l} indicate the
area fraction (quadrature weight) associated with node
*l*. **F**_{l} is
computed from **ϕ**_{l} using
a finite difference approximation to the fiber force density; see refs.
38 and 34 for details, including relevant finite
difference formulae.

The deformations, stresses, and resultant forces associated with
the aortic wall are approximated on a volumetric Lagrangian mesh. Let
*e* index the elements of this mesh, with
*V _{e}* indicating the volume associated
with element

$$\mathit{\chi}(\mathbf{X},t)=\sum _{m}{\mathit{\chi}}_{m}(t){\phi}_{m}(\mathbf{X})$$

(22)

$$\mathbf{G}(\mathbf{X},t)=\sum _{m}{\mathbf{G}}_{m}(t){\phi}_{m}(\mathbf{X}),$$

(23)

in which
**χ**_{m}(*t*)
are the nodal positions and
**G**_{m}(*t*) are
the nodal force densities. The deformation gradient is
computed by directly differentiating the approximation to
**χ**(**X**, *t*), and
is evaluated from the approximation to . The nodal
values **G**_{m}(*t*)
are determined so that

$${\int}_{V}\mathbf{G}(\mathbf{X},t){\phi}_{n}(\mathbf{X})\mathrm{d}\mathbf{X}=-{\int}_{V}\mathbb{P}(\mathbf{X},t){\nabla}_{\mathbf{X}{\phi}_{n}}(\mathbf{X})\mathrm{d}\mathbf{X}$$

(24)

for all *n*. This leads to a linear
system of equations that defines the nodal force densities in terms of
the nodal deformations. In practice, these integrals are evaluated
element-by-element using Gaussian quadrature rules that are exact for
the left-hand side but are generally only approximate for the right-hand
side.

To couple the Lagrangian and Eulerian discretizations, we employ
approximations to the integral transforms of the continuum equations
that replace the singular Dirac delta function δ(**x**)
= δ(*x*) δ(*y*)
δ(*z*) by a regularized delta function
δ_{h}(**x**) =
δ_{h}(*x*)
δ_{h}(*y*)
δ_{h}(*z*). In
our computations, we construct the three-dimensional regularized delta
function either by using the four-point one-dimensional regularized
delta function of Peskin [58], or
by using a broadened version of this function that has a spatial extent
of 8 meshwidths. The same regularized delta function is used for both
the thin model of the valve leaflets and the volumetric model of the
vessel wall.

For the leaflet model, we employ a Lagrangian-Eulerian coupling
approach frequently used with the IB method. The Lagrangian forces
**F** = (*F ^{x}, F^{y},
F^{z}*) associated with the leaflets are
converted into equivalent Eulerian forces

$${f}_{i-\frac{1}{2},j,k}^{x}=\sum _{l}{F}_{l}^{x}{\delta}_{h}({\mathbf{x}}_{i-\frac{1}{2},j,k}-{\mathit{\varphi}}_{l})\Delta {\mathbf{q}}_{l},$$

(25)

$${f}_{i,j-\frac{1}{2},k}^{y}=\sum _{l}{F}_{l}^{y}{\delta}_{h}({\mathbf{x}}_{i,j-\frac{1}{2},k}-{\mathit{\varphi}}_{l})\Delta {\mathbf{q}}_{l},$$

(26)

$${f}_{i,j,k-\frac{1}{2}}^{z}=\sum _{l}{F}_{l}^{z}{\delta}_{h}({\mathbf{x}}_{i,j,k-\frac{1}{2}}-{\mathit{\varphi}}_{l})\Delta {\mathbf{q}}_{l}.$$

(27)

This amounts to using the trapezoidal rule to discretize the integral transform (3). We employ the compact notation

$$\mathbf{f}={\mathit{\mathcal{S}}}^{\mathrm{l}}[\mathit{\varphi}]\mathbf{F}$$

(28)

to express this *force-spreading*
operation. The Eulerian velocity **u** = (*u*,
υ, *w*) is used to determine the dynamics of the
physical positions **ϕ** =
(ϕ^{x},
ϕ^{y},
ϕ^{z}) of the material points
of the valve leaflets via

$$\frac{\mathrm{d}}{\mathrm{d}t}{\varphi}_{l}^{x}=\sum _{i,j,k}{u}_{i-\frac{1}{2},j,k}{\delta}_{h}({\mathbf{x}}_{i-\frac{1}{2},j,k}-{\mathit{\varphi}}_{l})\Delta \mathbf{x},$$

(29)

$$\frac{\mathrm{d}}{\mathrm{d}t}{\varphi}_{l}^{y}=\sum _{i,j,k}{\upsilon}_{i,j-\frac{1}{2},k}{\delta}_{h}({\mathbf{x}}_{i,j-\frac{1}{2},k}-{\mathit{\varphi}}_{l})\Delta \mathbf{x},$$

(30)

$$\frac{\mathrm{d}}{\mathrm{d}t}{\varphi}_{l}^{z}=\sum _{i,j,k}{w}_{i,j,k-\frac{1}{2}}{\delta}_{h}({\mathbf{x}}_{i,j,k-\frac{1}{2}}-{\mathit{\varphi}}_{l})\Delta \mathbf{x}.$$

(31)

We use the notation

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathit{\varphi}={\mathit{\mathcal{R}}}^{\mathrm{l}}[\mathit{\varphi}]\mathbf{u}$$

(32)

to express this *velocity-restriction*
operation. Notice that the force-spreading and velocity-interpolation
operators are adjoints, i.e.,
^{l}[**ϕ**] =
^{l}[**ϕ**]*.

The Lagrangian-Eulerian coupling scheme used for the vessel
model is based on an approach developed by Griffith and Luo [39]. This approach is first to
construct a force-spreading operator, and then to determine a
velocity-restriction operator that is the adjoint of the force-spreading
operator. Briefly, for each element *e*, quadrature
points *Q*(*e*) are determined, and for
each quadrature point *q*
*Q*(*e*), ${\mathbf{X}}_{q}^{e}$ indicates the reference coordinates of the quadrature
point, and ${\omega}_{q}^{e}$ indicates the quadrature weight associated. The
Lagrangian force density **G** = (*G ^{x},
G^{y}, G^{z}*) is converted into the
equivalent Eulerian force density

$${g}_{i-\frac{1}{2},j,k}^{x}=\sum _{e}\sum _{q\in Q(e)}{G}^{x}({\mathbf{X}}_{q}^{e},t){\delta}_{h}({\mathbf{x}}_{i-\frac{1}{2},j,k}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t)){\omega}_{q}^{e},$$

(33)

$${g}_{i,j-\frac{1}{2},k}^{y}=\sum _{e}\sum _{q\in Q(e)}{G}^{y}({\mathbf{X}}_{q}^{e},t){\delta}_{h}({\mathbf{x}}_{i,j-\frac{1}{2},k}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t)){\omega}_{q}^{e},$$

(34)

$${g}_{i,j,k-\frac{1}{2}}^{z}=\sum _{e}\sum _{q\in Q(e)}{G}^{z}({\mathbf{X}}_{q}^{e},t){\delta}_{h}({\mathbf{x}}_{i,j,k-\frac{1}{2}}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t)){\omega}_{q}^{e}.$$

(35)

Notice that these formulae take advantage of the fact
that we may evaluate the approximations to **G** and
**χ** at arbitrary Lagrangian locations.
Specifically, we are not restricted to evaluating **G** and
**χ** at the nodes of the finite element mesh. In
practice, the quadrature rules are dynamically determined to ensure that
the physical spacing of the quadrature points is on average one-half of
a Cartesian meshwidth, even under the presence of extremely large
structural deformations. We use the notation

$$\mathbf{g}={\mathit{\mathcal{S}}}^{\mathrm{w}}[\mathit{\chi}]\mathbf{G}$$

(36)

to define the force-spreading operator. To determine the
corresponding velocity-restriction operator, we introduce a Lagrangian
velocity field **U** = (*U, V, W*), which is
defined in terms of nodal velocities
**U**_{m} via

$$\mathbf{U}(\mathbf{X},t)=\sum _{m}{\mathbf{U}}_{m}(t){\phi}_{m}(\mathbf{X}).$$

(37)

This velocity field is required to satisfy, for all
*n*,

$$\sum _{e}\sum _{q\in Q(e)}U({\mathbf{X}}_{q}^{e},t){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e}\phantom{\rule{0ex}{0ex}}=\sum _{e}\sum _{q\in Q(e)}(\sum _{i,j,k}{u}_{i-\frac{1}{2},j,k}{\delta}_{h}({\mathbf{x}}_{i-\frac{1}{2},j,k}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t))\Delta \mathbf{x}){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e},$$

(38)

$$\sum _{e}\sum _{q\in Q(e)}V({\mathbf{X}}_{q}^{e},t){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e}\phantom{\rule{0ex}{0ex}}=\sum _{e}\sum _{q\in Q(e)}(\sum _{i,j,k}{\upsilon}_{i,j-\frac{1}{2},k}{\delta}_{h}({\mathbf{x}}_{i,j-\frac{1}{2},k}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t))\Delta \mathbf{x}){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e},$$

(39)

$$\sum _{e}\sum _{q\in Q(e)}W({\mathbf{X}}_{q}^{e},t){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e}\phantom{\rule{0ex}{0ex}}=\sum _{e}\sum _{q\in Q(e)}(\sum _{i,j,k}{w}_{i,j,k-\frac{1}{2}}{\delta}_{h}({\mathbf{x}}_{i,j,k-\frac{1}{2}}-\mathit{\chi}({\mathbf{X}}_{q}^{e},t))\Delta \mathbf{x}){\phi}_{n}({\mathbf{X}}_{q}^{e}){\omega}_{q}^{e}.$$

(40)

This yields a system of linear equations to be solved
for the nodal values of **U**. Notice that eqs. (38)–(40) correspond to a
discrete component-wise approximation to eq. (12). We use the notation

$$\frac{\mathrm{d}}{\mathrm{d}t}\mathit{\chi}=\mathbf{U}={\mathit{\mathcal{R}}}^{\mathrm{w}}[\mathit{\chi}]\mathbf{u}.$$

(41)

Notice that
^{w}[**χ**] is a
nonlocal operator in the sense that it requires the solution of a system
of linear equations. Although not shown here, it is the case that
^{w}[**χ**] =
^{w}[**χ**]*. See
Griffith and Luo [39] for further
discussion.

For the models of the leaflets and for the vessel wall, the
Lagrangian-Eulerian coupling scheme is designed to conserve force and
torque [39,58]. Moreover, because adjoint operators are used
to *spread* force and to *interpolate*
velocity, the scheme conserves energy at least in its semi-discrete
(continuous time but discrete space) formulation.

The time stepping scheme is similar to that of Griffith and Luo
[39]. Let
Δ*t* indicate the (uniform) time step size,
and let [*t ^{n},
t*

The time stepping scheme proceeds as follows. First, predicted
intermediate approximations to the structural deformations
**ϕ** and **χ** at time ${t}^{n+\frac{1}{2}}$ are determined via

$$\frac{{\mathit{\varphi}}^{n+\frac{1}{2}}-{\mathit{\varphi}}^{n}}{\Delta t/2}={\mathit{\mathcal{R}}}^{\mathrm{l}}[{\mathit{\varphi}}^{n}]{\mathbf{u}}^{n},$$

(42)

$$\frac{{\mathit{\chi}}^{n+\frac{1}{2}}-{\mathit{\chi}}^{n}}{\Delta t/2}={\mathit{\mathcal{R}}}^{\mathrm{w}}[{\mathit{\chi}}^{n}]{\mathbf{u}}^{n}.$$

(43)

This is an application of the forward Euler time stepping scheme. Lagrangian force densities ${\mathbf{F}}^{n+\frac{1}{2}}$ and ${\mathbf{G}}^{n+\frac{1}{2}}$ are determined from the predicted structure configurations, and these Lagrangian force densities are spread to the Cartesian grid via

$${\mathbf{f}}^{n+\frac{1}{2}}={\mathit{\mathcal{S}}}^{\mathrm{l}}[{\mathit{\varphi}}^{n+\frac{1}{2}}]{\mathbf{F}}^{n+\frac{1}{2}},$$

(44)

$${\mathbf{g}}^{n+\frac{1}{2}}={\mathit{\mathcal{S}}}^{\mathrm{w}}[{\mathit{\chi}}^{n+\frac{1}{2}}]{\mathbf{G}}^{n+\frac{1}{2}}.$$

(45)

Next, we solve the discretized incompressible Navier-Stokes
equations for **u**^{n+1} and ${p}^{n+\frac{1}{2}}$ via a Crank-Nicolson-Adams-Bashforth scheme,

$$\rho (\frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^{n}}{\Delta t}+{[\mathbf{u}\xb7{\nabla}_{h}\mathbf{u}]}^{(n+\frac{1}{2})})=-{\nabla}_{h}{p}^{n+\frac{1}{2}}+\mu {\nabla}_{h}^{2}\left(\frac{{\mathbf{u}}^{n+1}+{\mathbf{u}}^{n}}{2}\right)+{\mathbf{f}}^{n+\frac{1}{2}}+{\mathbf{g}}^{n+\frac{1}{2}},$$

(46)

$${\nabla}_{h}\xb7{\mathbf{u}}^{n+1}=0,$$

(47)

with

$${[\mathbf{u}\xb7{\nabla}_{h}\mathbf{u}]}^{(n+\frac{1}{2})}=\frac{3}{2}{\mathbf{u}}^{n}\xb7{\nabla}_{h}{\mathbf{u}}^{n}-\frac{1}{2}{\mathbf{u}}^{n-1}\xb7{\nabla}_{h}{\mathbf{u}}^{n-1}.$$

(48)

In the uniform grid case, the discrete operators
_{h},
_{h} ·, and ${\nabla}_{h}^{2}$ are standard second-order accurate staggered-grid
finite difference approximations to the gradient, divergence, and
Laplace operators, respectively [33]. We compute **u** ·
_{h}**u** via a
staggered-grid version of the piecewise parabolic method (PPM) [33]. On locally-refined Cartesian
grids, these discrete operators are straightforward extensions of the
uniform grid operators [34].
Solving eqs. (46) and
(47) for
**u**^{n+1} and ${p}^{n+\frac{1}{2}}$ requires the solution of the time-dependent Stokes
equations. The Stokes equations are solved using an iterative Krylov
method with a multigrid preconditioner derived from the classical
projection method [9, 33].

Finally, we determine approximations to the structural deformations at the end of the time step via

$$\frac{{\mathit{\varphi}}^{n+1}-{\mathit{\varphi}}^{n}}{\Delta t}={\mathit{\mathcal{R}}}^{\mathrm{l}}[{\mathit{\varphi}}^{n+\frac{1}{2}}]{\mathbf{u}}^{n+\frac{1}{2}},$$

(49)

$$\frac{{\mathit{\chi}}^{n+1}-{\mathit{\chi}}^{n}}{\Delta t}={\mathit{\mathcal{R}}}^{\mathrm{w}}[{\mathit{\chi}}^{n+\frac{1}{2}}]{\mathbf{u}}^{n+\frac{1}{2}},$$

(50)

with

$${\mathbf{u}}^{n+\frac{1}{2}}=\frac{{\mathbf{u}}^{n+1}+{\mathbf{u}}^{n}}{2}.$$

(51)

This is an application of the explicit midpoint rule to the structural dynamics.

In the initial time step, we do not have a value for
**u**^{n−1} as required to
evaluate the Adams-Bashforth approximation to the convective term.
Consequently, in the initial time step, we employ a two-step Runge-Kutta
scheme.

Within each time step, coupling between the detailed fluid-structure interaction model and the reduced Windkessel model of the circulation is done as described previously [34,38]. In this approach, we perform only a single Stokes solve per time step, except in the initial time step, when we perform two Stokes solves within the context of a two-step Runge-Kutta scheme.

^{1}Keeping the loaded configuration fixed, if we assume a larger initial
pressure, then the result will be larger deformations from the computed unloaded
configuration to the prescribed loaded configuration, thereby leading to a
stiffer model at the prescribed loaded configuration. If we assume a smaller
initial pressure, then there will be smaller strains from the unloaded to the
loaded configuration, and thus a more compliant model at the target
configuration.

^{2}Similar convergence studies, although for substantially different
applications, are described in refs. 37
and 29.

Vittoria Flamini, Department of Mechanical and Aerospace Engineering, New York University Tandon School of Engineering, Brooklyn, New York, USA.

Abe DeAnda, Division of Cardiothoracic Surgery, Department of Surgery, University of Texas Medical Branch, Galveston, TX.

Boyce E. Griffith, Departments of Mathematics and Biomedical Engineering and McAllister Heart Institute, Phillips Hall, Campus Box 3250, University of North Carolina, Chapel Hill, North Carolina, USA, Phone: (919) 962-7110, Email: ude.cnu.liame@gecyob.

1. Alastrué V, Garía A, Peña E, Rodríguez J, Martínez M, Doblaré M. Numerical framework for patient-specific computational modelling
of vascular tissue. Int J Numer Method Biomed Eng. 2010;26(1):35–51.

2. Alastrué V, Peña E, Martínez MÁ, Doblaré M. Assessing the use of the opening angle method to enforce residual
stresses in patient-specific arteries. Ann Biomed Eng. 2007;35(10):1821–1837. [PubMed]

3. Azadani AN, Chitsaz S, Matthews PB, Jaussaud N, Leung J, Tsinman T, Ge L, Tseng EE. Comparison of mechanical properties of human ascending aorta and
aortic sinuses. Ann Thorac Surg. 2012;93(1):87–94. [PubMed]

4. Bellhouse BJ, Bellhouse FH. Mechanism of closure of the aortic valve. Nature. 1968;217:86–87. [PubMed]

5. Bols J, Degroote J, Trachet B, Verhegghe B, Segers P, Vierendeels J. A computational method to assess in vivo stresses and unloaded
configuration of patient-specific blood vessels. J Comput Appl Math. 2012

6. Bonow RO, Carabello BA, Chatterjee K, de Leon AC, Faxon DP, Freed MD, Gaasch WH, Lytle BW, Nishimura RA, O’Gara PT, O’Rourke RA, Otto CM, Shah PM, Shanewise JS. ACC/AHA 2006 guidelines for the management of patients with
valvular heart disease. Circulation. 2006;114(5):E84–E231. [PubMed]

7. Borazjani I. Fluidstructure interaction, immersed boundary-finite element
method simulations of bio-prosthetic heart valves. Comput Meth Appl Mech Eng. 2013;257:103–116.

8. Borazjani I, Ge L, Sotiropoulos F. Curvilinear immersed boundary method for simulating fluid
structure interaction with complex 3D rigid bodies. J Comput Phys. 2008;227(16):7587–7620. [PMC free article] [PubMed]

9. Cai M, Nonaka A, Bell J, Griffith B, Donev A. Efficient variable-coefficient finite-volume Stokes
solvers. Submitted.

10. Cardamone L, Valentín A, Eberth JF, Humphrey JD. Origin of axial prestretch and residual stress in
arteries. Biomech Model Mechanobiol. 2009;8(6):431–446. [PMC free article] [PubMed]

11. Carr JA, Savage EB. Aortic valve repair for aortic insufficiency in adults: A
contemporary review and comparison with replacement
techniques. Eur J Cardio Thorac Surg. 2005;25(1):6–15. [PubMed]

12. Cheng A, Dagum P, Miller DC. Aortic root dynamics and surgery: from craft to
science. Phil Trans R Soc B. 2007;362(1484):1407–1419. [PMC free article] [PubMed]

13. Cheng R, Lai YG, Chandran KB. Three-dimensional fluid-structure interaction simulation of
bileaflet mechanical heart valve flow dynamics. Ann Biomed Eng. 2004;32(11):1471–1483. [PMC free article] [PubMed]

14. Conti CA, Votta E, Della Corte A, Del Viscovo L, Bancone C, Cotrufo M, Redaelli A. Dynamic finite element analysis of the aortic root from
mri-derived parameters. Med Eng Phys. 2010;32(2):212–221. [PubMed]

15. Creane A, Kelly DJ, Lally C. Patient specific computational modeling in cardiovascular
mechanics. In: Lopez BC, Peña E, editors. Patient-Specific Computational Modeling. Springer; 2012. pp. 61–79.

16. Croft LR, Mofrad MRK. Computational modeling of aortic heart valves. In: De S, Guilak F, Mofrad MRK, editors. Computational Modeling in Biomechanics. Springer; 2010. pp. 221–252.

17. Crosetto P, Reymond P, Deparis S, Kontaxakis D, Stergiopulos N, Quarteroni A. Fluid–structure interaction simulation of aortic blood
flow. Comput Fluid. 2011;43(1):46–57.

18. Dagum P, Green GR, Nistal FJ, Daughteres GT, Timek TA, Foppiano LE, Bolger AF, Ingels NB, Miller DC. Deformational dynamics of the aortic root: modes and physiologic
determinants. Circulation. 1999;100(2):II54–II62. [PubMed]

19. De Hart J, Baaijens FPT, Peters GWM, Schreurs PJG. A computational fluid-structure interaction analysis of a
fiber-reinforced stentless aortic valve. J Biomech. 2003;36(5):699–712. [PubMed]

20. de Putter S, Wolters BJBM, Rutten MCM, Breeuwer M, Gerritsen FA, van de Vosse FN. Patient-specific initial wall stress in abdominal aortic
aneurysms with a backward incremental method. J Biomech. 2007;40(5):1081–1090. [PubMed]

21. Delfino A, Stergiopulos N, Moore JE, Meister JJ. Residual strain effects on the stress field in a thick wall
finite element model of the human carotid bifurcation. J Biomech. 1997;30(8):777–786. [PubMed]

22. Driscol TE, Eckstein RW. Systolic pressure gradients across the aortic valve and in the
ascending aorta. Am J Physiol. 1965;209(3):557–563. [PubMed]

23. Dumont K, Stijnen JMA, Vierendeels J, Van De Vosse FN, Verdonck PR. Validation of a fluid–structure interaction model of a
heart valve using the dynamic mesh method in fluent. Comput Meth Biomech Biomed Eng. 2004;7(3):139–146. [PubMed]

24. Esmaily Moghadam M, Bazilevs Y, Hsia TY, Vignon-Clementel IE, Marsden AL. A comparison of outlet boundary treatments for prevention of
backflow divergence with relevance to blood flow simulations. Comput Mech. 2011;48(3):277–291.

25. Fai TG, Griffith BE, Mori Y, Peskin CS. Immersed boundary method for variable viscosity and variable
density problems using fast linear solvers. I: Numerical method and
results. SIAM J Sci Comput. 2013;35(5):B1132–B1161.

26. Freeman RV, Otto CM. Spectrum of calcific aortic valve disease: Pathogenesis, disease
progression, and treatment strategies. Circulation. 2005;111(24):3316–3326. [PubMed]

27. Fung YC. What principle governs the stress distribution in living
organs? In: Fung YC, Fukada E, Wang J, editors. Biomechanics in China. Japan and USA: Science Press; 1983.

28. Fung YC. What are the residual stresses doing in our blood
vessels? Ann Biomed Eng. 1991;19(3):237–249. [PubMed]

29. Gao H, Wang HM, Berry C, Luo XY, Griffith BE. Quasi-static image-based immersed boundary-finite element model
of left ventricle under diastolic loading. Int J Numer Meth Biomed Eng. 2014;30(11):1199–1222. [PMC free article] [PubMed]

30. Gasser TC, Ogden RW, Holzapfel GA. Hyperelastic modelling of arterial layers with distributed
collagen fibre orientations. J R Soc Interface. 2006;3(6):15–35. [PMC free article] [PubMed]

31. Gee MW, Reeps C, Eckstein HH, Wall WA. Prestressing in finite deformation abdominal aortic aneurysm
simulation. J Biomech. 2009;42(11):1732–1739. [PubMed]

32. Govindjee S, Mihalic PA. Computational methods for inverse deformations in
quasi-incompressible finite elasticity. Int J Numer Method Eng. 1998;43(5):821–838.

33. Griffith BE. An accurate and efficient method for the incompressible
Navier-Stokes equations using the projection method as a
preconditioner. J Comput Phys. 2009;228(20):7565–7595.

34. Griffith BE. Immersed boundary model of aortic heart valve dynamics with
physiological driving and loading conditions. Int J Numer Method Biomed Eng. 2012;28(3):317–345. [PubMed]

35. Griffith BE, Hornung RD, McQueen DM, Peskin CS. An adaptive, formally second order accurate version of the
immersed boundary method. J Comput Phys. 2007;223(1):10–49.

36. Griffith BE, Hornung RD, McQueen DM, Peskin CS. Parallel and adaptive simulation of cardiac fluid
dynamics. In: Parashar M, Li X, editors. Advanced Computational Infrastructures for Parallel and Distributed
Adaptive Applications. Hoboken, NJ, USA: John Wiley and Sons; 2009.

37. Griffith BE, Lim S. Simulating an elastic ring with bend and twist by an adaptive
generalized immersed boundary method. Commun Comput Phys. 2012;12(2):433–461.

38. Griffith BE, Luo X, McQueen DM, Peskin CS. Simulating the fluid dynamics of natural and prosthetic heart
valves using the immersed boundary method. Int J Appl Mech. 2009;1(01):137–177.

39. Griffith BE, Luo XY. Hybrid finite difference/finite element version of the immersed
boundary method. (submitted)

40. Guy RD, Phillip B, Griffith BE. Geometric multigrid for an implicit-time immersed boundary
method. Adv Comput Math. 2015;41(3):635–662.

41. Humphrey JD. Cardiovascular Solid Mechanics: Cells, Tissues, and Organs. Springer Verlag; 2002.

42. Kim Y, Peskin CS. Penalty immersed boundary method for an elastic boundary with
mass. Phys Fluid. 2007;19:18 pages. 053,103.

43. Kim Y, Zhu L, Wang X, Peskin CS. On various techniques for computer simulation of boundaries with
mass. In: Bathe KJ, editor. Proceedings of the Second MIT Conference on Computational Fluid and
Solid Mechanics; Elsevier; 2003. pp. 1746–1750.

44. Kovács SJ, McQueen DM, Peskin CS. Modelling cardiac fluid dynamics and diastolic
function. Phil Trans R Soc Lond A. 2001;359(1783):1299–1314.

45. Lansac E, Lim HS, Shomura Y, Lim KH, Rice NT, Goetz W, Acar C, Duran CMG. A four-dimensional study of the aortic root
dynamics. Eur J Cardio Thorac Surg. 2002;22(4):497–503. [PubMed]

46. Lu J, Zhou X, Raghavan ML. Inverse elastostatic stress analysis in pre-deformed biological
structures: demonstration using abdominal aortic aneurysms. J Biomech. 2007;40(3):693–696. [PubMed]

47. Luo XY, Griffith BE, Ma XS, Yin M, Wang TJ, Liang CL, Watton PN, Bernacca GM. Effect of bending rigidity in a dynamic model of a polyurethane
prosthetic mitral valve. Biomech Model Mechanobiology. 2012;11(6):815–827. [PubMed]

48. Ma XS, Gao H, Griffith BE, Berry C, Luo XY. Image-based fluid-structure interaction model of the human mitral
valve. Comput Fluid. 2013;71:417–425.

49. Marom G, Haj-Ali R, Raanani E, Schäfers HJ, Rosenfeld M. A fluid–structure interaction model of the aortic valve
with coaptation and compliant aortic root. Med Biol Eng Comput. 2012;50(2):173–182. [PubMed]

50. May-Newman K, Lam C, Yin FCP. A hyperelastic constitutive law for aortic valve
tissue. J Biomech Eng. 2009;131(8):7 pages. 081,009. [PubMed]

51. McQueen DM, Peskin CS. Shared-memory parallel vector implementation of the immersed
boundary method for the computation of blood flow in the beating mammalian
heart. J Supercomput. 1997;11(3):213–236.

52. McQueen DM, Peskin CS. A three-dimensional computer model of the human heart for
studying cardiac fluid dynamics. Comput Graph. 2000;34(1):56–60.

53. McQueen DM, Peskin CS. Heart simulation by an immersed boundary method with formal
second-order accuracy and reduced numerical viscosity. In: Aref H, Phillips JW, editors. Mechanics for a New Millennium, Proceedings of the 20th International
Conference on Theoretical and Applied Mechanics (ICTAM 2000); Kluwer Academic Publishers.2001.

54. Mori Y, Peskin CS. Implicit second order immersed boundary methods with boundary
mass. Comput Meth Appl Mech Eng. 2008;197(25–28):2049–2067.

55. Murgo JP, Westerhof N, Giolma JP, Altobelli SA. Aortic input impedance in normal man: relationship to pressure
wave forms. Circulation. 1980;62(1):105–116. [PubMed]

56. Nichols WW, O’Rourke MF. McDonald’s Blood Flow in Arteries: Theoretical, Experimental,
and Clinical Principles. CRC Press; 2011.

57. Peskin CS. Flow patterns around heart valves: a numerical
method. J Comput Phys. 1972;10(2):252–271.

58. Peskin CS. The immersed boundary method. Acta Numer. 2002;11(0):479–517.

59. Peskin CS, McQueen DM. Mechanical equilibrium determines the fractal fiber architecture
of aortic heart valve leaflets. Am J Physiol Heart Circ Physiol. 1994;266(1):H319–H328. [PubMed]

60. Peskin CS, McQueen DM. Fluid dynamics of the heart and its valves. In: Othmer HG, Adler FR, Lewis MA, Dallon JC, editors. Case Studies in Mathematical Modeling: Ecology, Physiology, and Cell
Biology. Englewood Cliffs, NJ, USA: Prentice-Hall; 1996. pp. 309–337.

61. Reul H, Vahlbruch A, Giersiepen M, Schmitz-Rode TH, Hirtz V, Effert S. The geometry of the aortic root in health, at valve disease and
after valve replacement. J Biomech. 1990;23(2):181–191. [PubMed]

62. Roy S, Heltai L, Costanzo F. Benchmarking the immersed finite element method for
fluid-structure interaction problems ArXiv preprint
arXiv:1306.0936

63. Sacks MS, Yoganathan AP. Heart valve function: a biomechanical perspective. Philosophical Transactions of the Royal Society B: Biological
Sciences. 2007;362(1484):1369–1391. [PMC free article] [PubMed]

64. Sauren AAHJ. Ph.D. thesis. Technische Universiteit Eindhoven; 1981. The mechanical behaviour of the aortic valve.

65. Sellier M. An iterative method for the inverse elasto-static
problem. J Fluid Struct. 2011;27(8):1461–1470.

66. Shunk KA, Garot J, Atalar E, Lima JA. Transesophageal magnetic resonance imaging of the aortic arch and
descending thoracic aorta in patients with aortic
atherosclerosis. Journal of the American College of Cardiology. 2001;37(8):2031–2035. [PubMed]

67. Singh IM, Shishehbor MH, Christofferson RD, Tuzcu EM, Kapadia SR. Percutaneous treatment of aortic valve stenosis. Cleve Clin J Med. 2008;75(11):805–812. [PubMed]

68. Sotiropoulos F, Borazjani I. A review of state-of-the-art numerical methods for simulating
flow through mechanical heart valves. Med Biol Eng Comput. 2009;47(3):245–256. [PMC free article] [PubMed]

69. Stergiopulos N, Westerhof BE, Westerhof N. Total arterial inertance as the fourth element of the windkessel
model. Am J Physiol Heart Circ Physiol. 1999;276(1):H81–H88. [PubMed]

70. Swanson WM, Clark RE. Dimensions and geometric relationships of the human aortic valve
as a function of pressure. Circ Res. 1974;35(6):871–882. [PubMed]

71. Thubrikar M. The Aortic Valve. CRC Press; 1989.

72. Vaishnav RN, Vossoughi J. Estimation of residual strain in aortic segment. In: Hall CW, editor. Biomedical Engineering II: Recent Developments. Pergamon Press; 1983.

73. Vavourakis V, Papaharilaou Y, Ekaterinaris JA. Coupled fluid–structure interaction hemodynamics in a
zero-pressure state corrected arterial geometry. J Biomech. 2011;44(13):2453–2460. [PubMed]

74. Viscardi F, Vergara C, Antiga L, Merelli S, Veneziani A, Puppini G, Faggian G, Mazzucco A, Luciani GB. Comparative finite element model analysis of ascending aortic
flow in bicuspid and tricuspid aortic valve. Artif Organs. 2010;34(12):1114–1120. [PubMed]

75. Weinberg EJ, Kaazempur Mofrad MR. A multiscale computational comparison of the bicuspid and
tricuspid aortic valves in relation to calcific aortic
stenosis. J Biomech. 2008;41(16):3482–3487. [PubMed]

76. Westerhof N, Stergiopulos N, Noble MIM. Snapshots of Hemodynamics: An Aid for Clinical Research and Graduate
Education. Springer; 2010.

77. Wittek A, Karatolios K, Bihari P, Schmitz-Rixen T, Moosdorf R, Vogt S, Blase C. In vivo determination of elastic properties of the human aorta
based on 4d ultrasound data. Journal of the mechanical behavior of biomedical materials. 2013;27:167–183. [PubMed]

78. Yao J, Liu GR, Narmoneva DA, Hinton RB, Zhang ZQ. Immersed smoothed finite element method for
fluid–structure interaction simulation of aortic
valves. Comput Mech. 2012;50(6):789–804.

79. Zhu L, Peskin CS. Simulation of a flapping flexible filament in a flowing soap film
by the immersed boundary method. J Comput Phys. 2002;179(2):452–468.

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