|Home | About | Journals | Submit | Contact Us | Français|
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 . 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 . 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 . The compliance of the sinuses also helps to ensure the proper closure of the aortic valve . 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 . 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 .
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 . 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 . 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  and by Vaishnav and Vossoughi , 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. , 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 . This method was extended by Lu et al.  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.  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. , 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 . A more generic approach to including residual stress in arteries is described in Alastrué et al.  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  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.  and by Weinberg and Mofrad , 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 . Attempts to tackle this limitation have included restricting the problem to a two-dimensional analysis  and the implementation of complex remeshing algorithms in three-dimensional analyses .
An alternative approach to modeling interactions between the blood and the aortic valve and root is offered by Peskin’s immersed boundary (IB) method , which was introduced to simulate cardiac valve dynamics  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.  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. . 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  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  provides downstream loading for the aortic root, and a pressure waveform, taken from the clinical data set used to parameterize this reduced circulation model , 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:
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 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
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 C0 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
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 . 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 Es and a bending energy Eb,
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 cb is a spatially inhomogeneous bending stiffness and 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 cb model a thick, stenotic valve, whereas smaller values of cb 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 . 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 . 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  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.  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 . 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.  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.  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,
in which I1 = I1() = 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
in which I3 = I3() = 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.  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 xm denote positions of the finite element mesh nodes in the prescribed loaded configuration, let denote the corresponding nodal positions in the computed reference configuration after i iterations of the backward displacement algorithm, and let denote nodal positions in the computed loaded configuration when Xi is used as the unloaded reference configuration. Notice that xm, , and are all defined at the nodes of the finite element mesh.
The algorithm is straightforward: First, we initialize X0 := x (i.e., we use the deformed coordinates as an initial guess for the unloaded configuration). Then, given Xi, Xi+1 is determined by
Notice that if we define the displacement from the computed reference configuration to the deformed coordinates at step i by Di = χi − Xi, then (19) is equivalent to
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 L2 norm of the residual x − χi, i.e.
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.  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.  fit to the clinical data of Murgo et al. . 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 , a small amount of flow stabilization is applied in the fluid region within a distance of approximately 2h 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 , 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  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 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 ~ h4 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. , 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.
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. . 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.
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 , 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.5h 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 .
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. .
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  used to fit the Windkessel model  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 . 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.  and by Sellier . 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. , 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 . 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 ; 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  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.  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 .
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  as well as experimentally constrained models of the elasticity of the aortic valve leaflets . Further, it is now well known that the valve leaflets are in fact layered structures . 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 . 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 , except that here we also employ a finite element-based description of the aortic wall, using an approach introduced by Griffith and Luo .
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 indicate the position of the center of grid cell (i, j, k), in which h is the Cartesian grid meshwidth and Δx = h3 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 , υ is approximated at the locations , and w is approximated at . 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 Fl indicate the current position and Lagrangian force density of node l, respectively, and let Δql indicate the area fraction (quadrature weight) associated with node l. Fl 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 Ve indicating the volume associated with element e, so that V = eVe, let m index the mesh nodes, and let m (X) indicate the interpolatory Lagrangian finite element basis function associated with node m. The structural deformation and elastic force density are approximated in a standard manner via
in which χm(t) are the nodal positions and Gm(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 Gm(t) are determined so that
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 , 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 = (Fx, Fy, Fz) associated with the leaflets are converted into equivalent Eulerian forces f = (fx, fy, fz) via
This amounts to using the trapezoidal rule to discretize the integral transform (3). We employ the compact notation
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
We use the notation
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 . 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), indicates the reference coordinates of the quadrature point, and indicates the quadrature weight associated. The Lagrangian force density G = (Gx, Gy, Gz) is converted into the equivalent Eulerian force density g = (gx, gy, gz) by discretizing eq. (9), with δ(x) replaced by δh(x), using this quadrature rule, i.e.
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
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 Um via
This velocity field is required to satisfy, for all n,
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  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 . Let Δt indicate the (uniform) time step size, and let [tn, tn+1] = [nΔt, (n + 1)Δt] indicate the nth time interval. In this subsection, superscript indices always indicate the time step number. The state variables u, ϕ, and χ are defined at integer time steps. The pressure p, which is not a state variable of the system, is defined at half-integer time steps. Approximations to the state variables and to derived quantities such as the Lagrangian forces are also evaluated at half-integer time steps.
The time stepping scheme proceeds as follows. First, predicted intermediate approximations to the structural deformations ϕ and χ at time are determined via
This is an application of the forward Euler time stepping scheme. Lagrangian force densities and are determined from the predicted structure configurations, and these Lagrangian force densities are spread to the Cartesian grid via
Next, we solve the discretized incompressible Navier-Stokes equations for un+1 and via a Crank-Nicolson-Adams-Bashforth scheme,
In the uniform grid case, the discrete operators h, h ·, and are standard second-order accurate staggered-grid finite difference approximations to the gradient, divergence, and Laplace operators, respectively . We compute u · hu via a staggered-grid version of the piecewise parabolic method (PPM) . On locally-refined Cartesian grids, these discrete operators are straightforward extensions of the uniform grid operators . Solving eqs. (46) and (47) for un+1 and 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
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 un−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.
1Keeping 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.
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.