Quasi-static image-based immersed boundary-finite element model of left ventricle under diastolic loading

Finite stress and strain analyses of the heart provide insight into the biomechanics of myocardial function and dysfunction. Herein, we describe progress toward dynamic patient-specific models of the left ventricle using an immersed boundary (IB) method with a finite element (FE) structural mechanics model. We use a structure-based hyperelastic strain-energy function to describe the passive mechanics of the ventricular myocardium, a realistic anatomical geometry reconstructed from clinical magnetic resonance images of a healthy human heart, and a rule-based fiber architecture. Numerical predictions of this IB/FE model are compared with results obtained by a commercial FE solver. We demonstrate that the IB/FE model yields results that are in good agreement with those of the conventional FE model under diastolic loading conditions, and the predictions of the LV model using either numerical method are shown to be consistent with previous computational and experimental data. These results are among the first to analyze the stress and strain predictions of IB models of ventricular mechanics, and they serve both to verify the IB/FE simulation framework and to validate the IB/FE model. Moreover, this work represents an important step toward using such models for fully dynamic fluid–structure interaction simulations of the heart. © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons, Ltd.


INTRODUCTION
Cardiac diseases have been a major public health burden in industrialized countries for more than a century and are now the leading causes of morbidity and mortality worldwide. Despite intense effort, however, the biomechanics of heart disease and heart failure remain incompletely understood. It is well known that alterations in myocardial stress and strain distributions can have a significant impact on maladaptive processes such as hypertrophy [1,2], but it is not possible to measure directly intramural stress distributions in patients. Three-dimensional stress and strain distributions are readily provided by computational models of the heart, however, and such models therefore could assist 1200 H. GAO ET AL. in patient risk stratification and inform clinical decision making [3]. An improved understanding of ventricular biomechanics is also important for optimizing medical therapies and surgical procedures aimed at restoring normal heart function [4].
The development of realistic computational models of the heart is complicated by several factors. The heart is a multiphysics system in which the contractions of the myocardium are stimulated and coordinated by the electrophysiology of the heart. The passive and active mechanical properties of the myocardium are anisotropic and nonlinear, and the heart experiences large deformations resulting both from diastolic filling and also from systolic contraction. In addition, the heart is fundamentally a fluid-structure system, and a fully dynamic computational model of the heart must account for the substantial inertia of both the muscular heart wall and the blood. Finally, there are significant intersubject differences in cardiac anatomy and physiology that must be incorporated into fully subject-specific computational models. For a recent survey of approaches to patient-specific modeling of cardiac biomechanics, see ref. [5] .
Because of the complex anatomical geometry, large deformations, and nonlinear material response of the heart, most detailed models of cardiac mechanics have used nonlinear finite element (FE) methods. Nash and Hunter [6] developed an FE framework for large deformation heart simulation using the pole zero constitutive law. A subsequent study by Stevens et al. [7] suggested that the pole-zero law might need to be reformulated to achieve a better separation of material parameters associated with deviatoric stresses. Vetter and McCulloch [8] modeled the rabbit left ventricle (LV) in diastole using a three-dimensional FE method, and a good agreement with experimental measurement was found. Using a three-dimensional finite deformation model, Usyk et al. [9] studied the effects of the laminar fiber architecture on regional stress and strain.
Since the 1990's, various studies have demonstrated that the layered organization of the ventricular myofibrils results in an orthotropic passive material response [10,11]. Relatively few studies have used orthotropic material laws for LV simulations, however. Instead, the most earlier work has described the myocardium as transversely isotropic [12,13]. The Holzapfel-Ogden model of the ventricular myocardium [10] provides a description of the orthotropic and hyperelastic behavior of cardiac muscle that is based on strain invariants that reflect the structure of the myocardium. Among the first organ-scale simulations to use this model were by Göktepe et al. [14], who used it with an FE model of the heart with an idealized biventricular geometry.
In earlier work, we also used the Holzapfel-Ogden strain-energy functional to model ventricular mechanics. We developed a realistic human left ventricular geometry derived from clinical magnetic resonance (MR) imaging data along with a rule-based fiber architecture. Using this model along with a nonlinear FE method for the structural mechanics, we studied the effects of changes in fiber organization and material properties on diastolic mechanics, and we demonstrated that changes to the fiber orientations have substantially larger effects on the passive response of the myocardium than changes to the sheet orientations [15]. We also showed that the effects of the residual stress, determined either from measured residual strains or by the opening-angle method, can be easily incorporated in a modified Holzapfel-Ogden model [16].
Dynamic fluid-structure interaction (FSI) models of the heart have also been developed. Perhaps the most widely used approach to construct such models is to use body-conforming discretizations of the fluid and structure, for example, via an arbitrary Lagrangian-Eulerian (ALE) formulation [17][18][19][20]. Watanabe et al. [17] developed an FSI model of the LV using an ALE FE method with strong coupling that also included a description of the electrophysiology of the LV. More recently, Nordsletten et al. [18] developed an FE framework for simulating blood flow and myocardial dynamics in both the diastolic and systolic phases of the cardiac cycle via a nonconforming ALE FE scheme. Their method has been used to study the effects of left ventricular assist devices on the heart [19] and to investigate the diastolic function of patients with hypoplastic left heart syndrome [20].
Fluid-structure interaction methods involving body-fitted grids generally require dynamic mesh generation for problems that involve large structural deformations and this requirement substantially complicates the implementation of such schemes. The immersed boundary (IB) method [21] offers an alternative to body-fitted methods that avoid the need for dynamic mesh generation. The IB method is a general approach to modeling systems in which an incompressible elastic structure is immersed in a viscous incompressible fluid. Within the IB framework, the momentum, viscosity, LV MECHANICS BY IB/FE: APPLICATION TO DIASTOLE 1201 and incompressibility of the coupled fluid-structure system are described in Eulerian form, and the structural deformations, stresses, and forces are described in Lagrangian form. Integral equations with Dirac delta function kernels couple the Eulerian and Lagrangian variables, and when the continuous equations are discretized, the singular delta function is replaced by a carefully constructed regularized version of the delta function. A key advantage of the IB approach to FSI is that it does not require the construction of conforming discretizations of the fluid and solid. This approach is therefore extremely well-suited for applications such as cardiac dynamics that involve large structural deformations.
Conventionally, many of the structural models used with the IB method have described the elasticity of the immersed body using systems of elastic fibers. Such fiber models have been used to develop a variety of biofluid dynamics models, including detailed three-dimensional models of the heart [22][23][24][25][26][27][28] and its valves [29][30][31][32]. Although fiber-based elasticity models are well-suited to describing anisotropic tissues such as the ventricular myocardium, it is challenging to use fiber models along with constitutive models, such as the Holzapfel-Ogden model [10], that account for the experimentally observed orthotropic properties of the ventricular myocardium. Further, in ventricular mechanics, stress and strain are important biomarkers [2], and despite extensive work on simulating cardiac dynamics by the IB method, we are aware of no previous studies that have assessed the accuracy of the stress and strain predictions of IB models of cardiac mechanics.
The IB method is not restricted to fiber-based elasticity models, and various extensions of the IB method have been introduced that use FE-based structural models [33][34][35][36][37][38][39]. In this study, we apply one of these schemes, a hybrid finite difference-finite element IB method [38] to simulate left ventricular biomechanics under diastolic loading conditions, and we present an initial verification of the accuracy of the predicted stress and strain distributions by comparing the model results to those generated by a well-established FE solver. As in our earlier work using a purely structural nonlinear FE method [15], here, we use the structure-based hyperelastic ventricular model of Holzapfel and Ogden [10] along with an anatomical geometry derived from clinical MR imaging data obtained from a healthy human subject and a rule-based myocardial fiber structure. Numerical predictions of the IB/FE LV model are shown to be in good qualitative and quantitative agreement with results from the earlier FE-based LV model, and the predictions of both models are shown to be in good agreement with previous computational and experimental data. We also demonstrate that improved numerical predictions of the intramural strain and stress can be obtained by employing straightforward modifications of the constitutive model that effectively redistribute part of the hydrostatic pressure from the Eulerian grid to the Lagrangian mesh.
Although the present study considers only the quasi-static mechanics of the LV under diastolic loading conditions, it is nonetheless among the first to employ the IB method to simulate cardiac mechanics along with realistic descriptions of cardiac anatomy and physiology. Moreover, our results serve to verify the IB/FE methodology and also provide an initial validation of this IB/FE model of ventricular mechanics. Our results also are among the first to demonstrate that the IB method is capable of producing quantitatively accurate stress and strain distributions for realistic models of the heart and is thereby a suitable framework for developing realistic FSI models of the heart.

METHODOLOGY
In this section, we give the details of developing an IB/FE LV model based on a subject-specific LV geometry from clinical imaging. We begin with a brief introduction of the IB formulation and the constitutive model of the passive myocardium. We then discuss modifications of the structural stress tensor that improve the accuracy of the computed structural stresses, as demonstrated empirically in Section 3. Next, we give an overview of the numerical methods, including the spatial discretization and time-stepping schemes. We subsequently describe the reconstruction of a human left ventricular model from MRI data. We conclude by describing the benchmark nonlinear FE LV model, which is implemented using the commercial ABAQUS nonlinear FE software and by providing a listing of the boundary and loading conditions of the IB/FE model and its numerical discretization parameters.

Immersed boundary formulation
In the IB approach to modeling FSI, the momentum, velocity, and incompressibility of the coupled fluid-structure system are described in Eulerian form, whereas the structural deformation, stress, and forces of the immersed solid body are described in Lagrangian form. In this work, the physical domain occupied by the fluid-structure system is denoted R 3 , in which x D .x 1 ; x 2 ; x 3 / 2 are fixed Eulerian (physical) coordinates, and the reference configuration of the immersed solid is denoted U R 3 , in which X D .X 1 ; X 2 ; X 3 / 2 U are fixed Lagrangian (material) coordinates. W .U; t / 7 ! is a time-dependent deformation mapping relating material and physical coordinates, so that .X; t/ 2 gives the physical position of material point X at time t . The physical region occupied by the immersed structure at time t is .U; t / D s .t / , and the physical region occupied by the fluid at time t is f .t / D n s .t /. The fluid-solid interface in the physical configuration at time t is @ s .t / and has an outward unit normal n.x; t/ (i.e., n is oriented pointing out the solid region s .t / into the fluid region f .t /). The boundary of the solid body in the fixed reference coordinate system is @U and has an outward unit normal N.X/.
The continuous equations of motion for the coupled fluid-structure system are [36] as follows: f.x; t/ D Z U r X P e .X; t/ ı.x .X; t// dX Z @U P e .X; t/ N.X/ ı.x .X; t// dA.X/; (4) in which u.x; t/ is the Eulerian velocity field that provides a spatial description of the velocity of the coupled fluid-structure system, p.x; t/ is the Eulerian pressure field that enforces the incompressibility of the coupled fluid-structure system, f.x; t/ is the Eulerian elastic force density, P e .X; t/ is the first Piola-Kirchhoff elastic stress tensor associated with the immersed solid, and ı.x/ D ı.x 1 / ı.x 2 / ı.x 3 / is the three-dimensional Dirac delta function. Implicit in the formulation (1)-(4) are the assumptions that the mass density and viscosity of the fluid and solid are equal. These assumptions are not essential, however, and versions of the IB method have been developed that permit the use of spatially varying structural mass densities [39][40][41][42][43][44] and viscosities [39,44]. The IB formulation includes two interaction equations that relate Lagrangian and Eulerian quantities via integral transforms with Dirac delta function kernels. The first interaction equation (3) determines the deformations of the immersed solid from the Eulerian velocity field u.x; t/. The presence of viscosity ensures that the Eulerian velocity field is continuous, and so Equation (3) implies @ @t .X; t/ D u. .X; t/; t/: This condition holds along the fluid-solid interface and therefore is equivalent to the velocity continuity (i.e., no slip) condition typically imposed at fluid-solid interfaces by ALE methods for FSI. Because r u.x; t/ D 0, the immersed solid is automatically treated as an exactly incompressible material. In particular, in the continuous equations, it is not necessary to impose the incompressibility constraint in the Lagrangian equations (e.g., via penalty terms in the elastic energy associated with the immersed solid). However, in practice, we empirically find that such penalty terms improve the stress predictions of the discrete LV model used in this work; see Section 3. The second interaction equation (4) converts the Lagrangian description of the elastic stresses generated by deformations of the immersed solid into a corresponding Eulerian force density. It can be shown [36] that for f defined by Equation (4), the right-hand side of the momentum equation (1) is equivalent to in which .x; t/ is the total Cauchy stress tensor of the coupled fluid-structure system, f .x; t/ is the stress tensor of a viscous incompressible fluid, and e .x; t/ is the stress tensor that corresponds to the elastic response of the immersed structure. The fluid-like stress tensor f is defined by in which I is the identity tensor, and the elastic stress tensor e is defined in terms of P e by e D´1 J P e F T for in which F D @ =@X is the deformation gradient tensor associated with the mapping W .U; t / 7 ! and J D det.F / is the Jacobian determinant. The stress tensors f and e are both defined for all x 2 , and although e is nonzero only in the solid region s .t / , fluid-like stresses are present throughout the physical domain. It is possible to remove the viscous stresses within the immersed body via an extension of the present formulation [39,44]. Such viscous stresses are appropriate to include in models of many biological tissues, including the ventricular myocardium, in which water comprises a substantial fraction of the tissue volume.
Because e is nonzero only within the solid region s .t /, continuity of traction at the fluid-solid interface generally requires discontinuities in the pressure and viscous stress along that interface. Such discontinuities are automatically accounted for in the IB formulation. Let indicate the value of a discontinuity in the bracketed quantity across an interface, which we define for a Eulerian function g.x/ at a position x 2 @ s .t / by in which g C .x/ and g .x/ indicate the limiting values of g.x 0 / as x 0 D x˙sn approaches x from f .t / and s .t /, respectively ‡ . Continuity of traction along @ s .t / requires n D f n C e n D 0; but the definition of e in Equation (8) implies that along the fluid-solid interface @ s .t /, e n D e n; (11) because e Á 0 in f .t /. Thus, along @ s .t /, there must be a corresponding discontinuity in the fluid-like traction, f n D e n: It can be shown that this implies p D n e n; for any unit vector t tangent to the interface. Thus, the limiting value of the elastic traction on @ s .t / determines the magnitudes of the discontinuities in the pressure and viscous stress. By Equation (8) and Nanson's relation, the limiting value of the elastic traction is given by in which T 1 and T 2 are mutually orthogonal unit tangent vectors to @U in the material coordinate system, and kF T 1 F T 2 k 1 is a conversion factor that relates densities per unit area in the reference and physical configurations at time t . P e N is precisely the force density that appears in the second integral transform in Equation (4), and in the IB formulation, this term imposes the required discontinuities in the pressure and viscous stress [45]. Moreover, these jump conditions are equivalent to the traction continuity conditions typically imposed at fluid-solid interfaces by ALE methods for FSI.

Orthotropic constitutive model
We describe the orthotropic elastic properties of the myocardium using a hyperelastic strain-energy functional introduced by Holzapfel and Ogden [10] W D a 2b expOEb.
in which a, b, a i , and b i (for i D f, s, and fs) are eight non-negative material parameters. Here, I 1 D tr.C/ is the first invariant of the right Cauchy-Green deformation tensor C D F T F , and I ? 4f , I ? 4s , and I 8fs are structure-based invariants that account for the anisotropy and shear properties of the myocardium. Specifically, let f 0 D f 0 .X/ and s 0 D s 0 .X/ denote the fiber and sheet axes in the reference configuration. In terms of these material direction fields, the basic structure-based invariants are defined by in which f D F f 0 and s D F s 0 are the deformed fiber and sheet axes in the physical configuration at time t . Thus, I 4f is the square of the stretch in the fiber direction, I 4s is the square of the stretch in the sheet direction, and I 8fs measures the relative orientation and deformation in the fiber and sheet directions. The constitutive model is not defined in terms of I 4f and I 4s but rather in terms of modified fiber and sheet invariants I ? 4i , which are given by for i D f and s. This ensures that the elastic energies and stresses associated with the collagen fibers are nonzero only in states of extension and not compression. For further discussion on this constitutive model, the reader is referred to the work of Holzapfel and Ogden [10]. The first Piola-Kirchhoff stress tensor is computed from W via P e D @W @F : In this work, we use the constitutive model parameters determined in our previous study [15],

Modified stress tensor
As previously discussed, nonzero values of P e N along the fluid-solid interface @U generate discontinuities in the pressure p and shear stress OEru C .ru/ T . These discontinuities are not artifacts of the IB formulation and are in fact required to obtain traction continuity at the fluid-solid interface. In the spatially discrete IB equations, such pressure discontinuities lead to spurious volume changes, and although these spurious fluxes vanish under grid refinement, they do so only at a firstorder rate [47]. It is straightforward to verify that with the present constitutive model, P e N ¤ 0 along the fluid-solid interface for F ¤ 0.
To reduce the magnitude of the pressure discontinuity in the Eulerian pressure field and thereby reduce the magnitude of spurious volume loss at fluid-solid interfaces, we have found that it is useful to employ a modified stress tensor Q P e defined so that Q P e D 0 for F D I. Specifically, we use so that in s .t /, The modified stress tensor Q P e differs from P e only by a modification to the pressure. Specifically, in the solid region s .t /, the physical pressure is in which This modified stress tensor is similar to modifications proposed for simpler neo-Hookean models to ensure that the stress vanishes in the reference configuration [48]. We briefly derive this modified stress tensor from a modified strain-energy functional, as follows. Let For an incompressible material, J D 1 and log.J / D 0, and therefore, Q W is identical to W in Equation (16). In addition, the definitions of W and Q W differ only in the isotropic term. Thus, it suffices to consider only the terms We compute the corresponding stresses as For J D 1, which is exactly imposed in the continuum setting but only approximately imposed in the discrete setting, we have This is precisely the form of the stress modification employed in this work. Although incompressibility is enforced in the Eulerian equations (1) and (2), when the continuous equations are discretized, the numerical interpolation of the Eulerian velocity to the solid region may not always yield a divergence-free discrete Lagrangian velocity field. This numerical error vanishes under grid refinement [47], but it can nonetheless lead to non-negligible errors in Lagrangian volume conservation at practical grid spacings, especially in three spatial dimensions. In the present context, such volume conservation errors appear to lead primarily to errors in the computed stress distributions rather than errors in the overall structural displacement. We obtain improved accuracy in the computed stress by including a penalty term to impose approximately the incompressibility constraint in Lagrangian form, namely, in which I 3 D det.C/ D J 2 andˇs > 0 is a penalty parameter. When using this modified stress, the physical pressure in the solid region s .t / is The Eulerian pressure field p is not a state variable of the system and automatically adjusts to account for the inclusion of additional pressure-like terms in the structural stress. Thus, at least in the continuum setting, the dynamics are identical regardless of whether we use P e or Q Q P e to determine the elastic stress. Except where otherwise noted, we use Q Q P e withˇs D5.0e6 dyne/cm 2 to determine the stress of the body, and we show that in the spatially discrete setting, the modified stress yields improved numerical accuracy (although not order of accuracy).

Numerical methods
The IB formulation (1)-(4) is not well-suited for discretization methods that employ standard FE methods for the structural equations. Consequently, in practice, we use a weak form of the Lagrangian equations: in which F.X; t/ is the Lagrangian elastic force density and V.X/ is an arbitrary Lagrangian test function that is not assumed to vanish on @U . This formulation permits us to approximate the Eulerian equations on a Cartesian grid and to approximate the Lagrangian equations on an unstructured hexahedral FE mesh that corresponds to the imaged LV geometry. We summarize our discretization approach in the following subsections.
2.4.1. Eulerian spatial discretization. Let .i; j; k/ label the grid cells of the Cartesian grid used to discretize the Eulerian equations. We denote the physical location of the center of cell .i; j; k/ by x i;j;k D ..i C 1=2/h; .j C 1=2/h; .k C 1=2/h/, in which h D x 1 D x 2 D x 3 is the Cartesian grid spacing. We use a staggered-grid discretization of the Eulerian velocity and force fields, so that each vector component is approximated at the center of the Cartesian cell face to which that component is normal. Specifically, the x 1 components of u and f are approximated at spatial locations x i 1 2 ;j;k D .ih; .j C 1=2/h; .k C 1=2/h/, the x 2 components are approximated at spatial locations x i;j 1 2 ;k D ..i C1=2/h; j h; .k C1=2/h/, and the x 3 components are approximated at spatial locations x i;j;k 1 2 D ..i C 1=2/h; .j C 1=2/h; kh/. The pressure p is approximated at the centers of the grid cells, that is, at spatial locations x i;j;k . Standard second-order accurate finite difference approximations to the divergence, gradient, and Laplace operators are employed [30,49]. The convective term is discretized using a version of the piecewise parabolic method [49][50][51].

Lagrangian spatial discretization.
Let l label the nodes of the FE mesh used to discretize the Lagrangian equations. The Lagrangian deformation mapping is approximated at the nodes of the FE mesh and is interpolated to arbitrary material coordinates using the FE shape functions, that is, in which l is the physical position of FE mesh node l and ' l is the nodal basis function associated with node l. The deformation gradient is approximated in the element interior by directly differentiating the approximation to in (37), and the stress P e is evaluated directly from the approximation to F . From P e , we compute an approximation to the Lagrangian body force F that satisfies for all Lagrangian basis functions ' m . This implicitly defines a system of linear equations for the nodal force densities F l involving the FE mass matrix. We use a standard selectively reduced integration approach to avoid volumetric locking [52], whereby reduced-order quadrature rules are used for the pressure normalization (p s ¤ 0) and penalty (ˇs ¤ 0) terms in the modified stress, whereas full-order quadrature is used for the remaining terms of the structural stress tensor.

Lagrangian-Eulerian interaction.
Our approach to discretize the interaction equations (34) and (35) is first to define a force-spreading operator S that determines f from F and then to define the velocity-restriction operator R that determines @ =@t from u to be the adjoint of the spreading operator, so that R D S . This construction yields a numerical scheme that, at least in the semi-discrete case, conserves energy during Lagrangian-Eulerian interaction [21]. In these approximations, we replace the singular delta function kernel with a regularized delta function of the tensor product form ı h .x/ D ı h .x 1 / ı h .x 2 / ı h .x 2 /. In this work, we construct the three-dimensional regularized delta function in terms of the four-point one-dimensional regularized delta function of Peskin [21].

H. GAO ET AL.
To convert the Lagrangian elastic force density into the corresponding Eulerian elastic force density, we compute f D .f 1 ; . .
in which F is defined by (38), and the right-hand side integrals are approximated via Gaussian quadrature. We use the notation f D SOE F; (43) in which the discrete operator S D SOE is defined by Equations (40)- (42).
To determine the velocity of the immersed structure, we first note that from (37), U D @ =@t satisfies that is, U is in the space of Lagrangian functions spanned by the nodal basis functions and has nodal values U l D d l =dt . It can be shown [38] that the operator R that determines U from u which satisfies R D S is implicitly defined by requiring U to satisfy Z for all basis functions ' m , in which the right-hand side of (45) is computed via the same quadrature rules as used in (38) and in which Notice that U IB is a continuous Lagrangian velocity field, but that U IB is not in the span of the nodal basis functions. From (44) and (45), it can be seen that U.X; t/ is the L 2 projection of U IB onto the space of functions spanned by the Lagrangian basis functions. We use the notation

Time stepping.
We use a time discretization that combines the explicit midpoint rule for the Lagrangian equations with a Crank Nicolson-Adams Bashforth scheme for the Eulerian equations. Let u n and n indicate the approximations to the Eulerian velocity field and Lagrangian deformation field at time t n D nt , in which t is the time step size, and let p nC 1 2 indicate the approximation to the Eulerian pressure field at time t nC 1 2 D .n C 1=2/t . The time discretization proceeds as follows.
First, we compute an approximation to at t nC 1 2 via We compute a corresponding approximation to the Lagrangian force density F nC 1 2 from nC 1 2 , and we spread this force to the Cartesian grid via Next, we solve the incompressible Navier-Stokes equations for u nC1 and in which OEu r h u nC 1 2 D 3 2 u n r h u n 1 2 u n 1 r h u n 1 : Solving (52) and (53) requires the solution of the incompressible Stokes equations, and to do so, we employ an iterative Krylov method preconditioned by a block-multigrid preconditioner based on a pressure-free projection method [49]. Finally, we determine the structure configuration at time t nC1 via nC1 n t D ROE nC 1 2 u nC1 C u n 2 : In the initial time step, we do not have lagged velocity values and therefore cannot compute the multistep approximation to the convective term used in subsequent time steps. Consequently, for the first time step, we instead use an explicit midpoint rule to approximate the convective term, which requires an additional Stokes solve. [27,28,53], which provides an adaptive and distributed-memory parallel infrastructure for developing FSI models that use IB method. IBAMR leverages functionality provided by other freely available software libraries, including SAMRAI [54,55], PETSc [56][57][58], and libMesh [59,60].

Imaging-derived human left ventricle model
A human left ventricular model is reconstructed from a cardiac MR imaging study on a healthy volunteer (male, age 28 years), performed at the British Heart Foundation MRI Facility at the University of Glasgow, as described previously [15]; see Figure 1 hexahedral LV mesh composed of 48,050 elements and 53,545 nodes; see Figure 1(b). A rule-based myocardial fiber generation procedure based on the work of Potse et al. [61] is employed to determine the fiber and sheet directions within the myocardium. In the present work, the fiber angle rotates from 60 ı to 60 ı from endocardium to epicardium, and the sheet angle rotates from 45 ı to 45 ı , thereby corresponding to a normal, healthy LV [11]; see Figure 1(c,d).

Structural analysis by ABAQUS
To verify the implementation of the IB/FE model, the same structural model is also analyzed using commercial ABAQUS (SIMULIA, Providence, RI, USA) nonlinear FE software, as in earlier work [15,16].

Physical and numerical parameters, boundary conditions, and loading conditions
In our simulations, we set D1.0 g/ml and D0.04 cP § , and is taken to be a 15 cm 15 cm 15 cm box that is discretized on a regular Cartesian grid. The integrals of the interaction equations are discretized using dynamically generated Gaussian quadrature rules that ensure a density of at least two quadrature points per Cartesian mesh width. A Cartesian grid convergence analysis is performed to demonstrate the sensitivity of the results on the Eulerian spatial mesh width. We consider N N N Cartesian grids for N D 64, 80, 96, 112, and 128. Unless otherwise noted, in the majority of our simulations, we use N D 96, which we find to yield essentially grid-converged results; see Section 3. In this case, the Eulerian grid spacing is h D x 1 D x 2 D x 3 D § In this study, we consider only the equilibrium configuration of the LV model under steady loading conditions, and although the values of and affect the transient dynamics, they have no effect on the steady-state results.

RESULTS
We first perform a grid sensitivity analysis for N N N Cartesian grids for N D 64, 80, 96, 112, and 128. Table I shows the differences in the computed displacements, fiber strains, and fiber stresses when compared with the ABAQUS-based model. Notice that only small differences in the computed values are observed for N > 96. Consequently, we use N D 96 for all subsequent IB/FE computations.
The displacement magnitudes determined by both models at a loading pressure of 7.5 mmHg are shown in Figure 3 panels (a) and (b), whereas panels (c) and (d) show the normed difference in the computed displacement. The maximum displacement obtained using the IB/FE model is 9.31 mm, which is within 1% of the maximum displacement of 9.40 mm obtained from the ABAQUS-based FE model. In addition, the spatial distributions of displacement are nearly identical in both models; compare panels (a) and (b) of Figure 3. Figure 3(c) plots the spatial distribution of the norm of the difference of the displacement vectors. In most of the volume, the differences are seen to be negligible, much less than 0.1 mm, that is, much less than a 1% difference. Figure 3(d) shows Table I. Cartesian grid sensitivity analysis, comparing the mean displacement, fiber strain, and fiber stress differences related to the results from the ABAQUS-based model. Notice that essentially grid-resolved results are obtained for N > 96.     Figure 4(a,b) for both the IB/FE and FE models. Figure 4(c,d) show the logarithmic fiber strain distributions for both the IB/FE and FE models. Stress and strain distributions generated by the two models are similar, although some minor differences can be seen. For instance, along the basal plane, there are modest differences, which likely result from differences in the treatment of displacement-type boundary conditions in the two models.
To facilitate a closer inspection of the results, seven selected transmural paths across the LV free wall are constructed, as shown in Figure 5. Figure 5(a-g) show fiber strain comparisons along the seven paths, and Figure 5(h) shows the differences between the IB/FE and FE-based models along those selected paths. The results indicate that the general fiber strain distributions obtained from both models are similar; however, the IB/FE model tends to yield slightly larger fiber strains close to the epicardium. The relative strain difference is generally less than 8%; see Figure 5(h). Figure 6(a-g) show fiber stress distributions along the same seven paths. Again, similar results are found in Figure 5: the fiber stress is similar in the two models, and the relative fiber stress difference is less than 8%. Figure 7 compares the results of the IB/FE-based and FE-based models to published strain distributions from earlier experimental and computational studies of the canine heart [12,63,64]. These data show that there is an overall qualitative and quantitative agreement between the IB/FE model predictions, the FE-based model, and the earlier experimental and computational results, especially for E rr (radial strain), E cr (circumferential-radial shear), and E lr (longitudinal-radial shear). The IB/FE-based and FE-based models are also in good quantitative agreement for E cc (circumferential strain), E ll (longitudinal strain), and E cl (circumferential-longitudinal shear), although in those cases, there are some discrepancies between our models and the results from earlier experimental and computational models. We remark that the strain distributions determined by the FE model in this work are moderately different from those determined in our earlier work [15], in which differ-1214 H. GAO ET AL. ent boundary conditions were imposed on the basal plane of the model LV. This suggests that these strain measurements are relatively sensitive to changes in boundary conditions. Figure 8 compares the end-diastolic pressure-volume relationship (EDPVR) of the IB/FE and FE models to EDPVR measurements from ex vivo human hearts by Klozt et al. [65]. These curves are normalized by EDV n D .EDV V 0 /=.V 30 V 0 /, in which V 0 is the initial unloaded LV volume, and V 30 is the LV volume with a loading of 30 mmHg. It is clear that there is good agreement between the IB/FE and FE model results and the data of Klozt et al. Figure 9 shows the effects of using the modified structural stress tensors, which include an additional pressure-like term p s , an incompressibility penalty term that is activated forˇs > 0, or both, on the average logarithmic fiber strain and the stress along the selected seven paths, as defined in Figure 5(i) and their differences. For the fiber strain, the differences are mainly concentrated at the endocardial and epicardial surfaces, where the IB/FE method yields lower-order accuracy because of its treatment of the Eulerian pressure discontinuities at those surfaces. A modified structural stress tensor obtained when using only the penalty term (i.e., p s D 0 andˇs ¤ 0) predicts similar fiber strains as when we include both pressure normalization and the penalty term (i.e., p s ¤ 0 anď modifications to structural stress tensors on the computed incompressibility of the myocardium. By includingˇs ¤ 0 in the structural stress tensor, J is very close to 1, indicating near incompressibility, whereas withˇs D 0, some elements can be compressed substantially. Table II suggests that for the discretized LV model used in this study, it is necessary to include terms that act to impose incompressibility in a Lagrangian sense even though the computed Eulerian velocity field is discretely divergence free. Figure 10 shows the velocity fields at t D 0:75 s and t D 5:0 s. At the beginning of the simulation, when the ventricular wall is initially pressurized, there is inflow into the LV as the chamber expands. By time t D 5:0 s, the model has essentially reached equilibrium and there is no further flow. Because of the lack of valves and simplified boundary conditions, the flow patterns from the present model are not physiological.

DISCUSSION
The IB method has been widely used to model cardiac FSI and to predict flow patterns and structural deformations; however, the fiber-based structural models conventionally used in IB simulations of cardiac dynamics are restricted to somewhat simple constitutive models. This makes it difficult to incorporate experimentally based material models into IB computations and thereby limits the ability of these models to predict tissue stress and strain distributions. Quantities such as the enddiastolic left-ventricular stress and strain distributions in patients with diastolic heart failure are  [64] and May-Newman et al. [63] and numerical results from a canine left ventricular model of Guccione et al. [12]. E cc , circumferential strain; E rr , radial strain; E ll , longitudinal strain; E cr , circumferential-radial shear; E cl , circumferential-longitudinal shear; E lr , longitudinal-radial shear.
important biomarkers and, moreover, are not accessible via standard clinical imaging. A typical approach to obtain these stress and strain data is to employ a nonlinear FE method along with a detailed model of the myocardium [10,[14][15][16]. Herein, we have demonstrated that the IB/FE framework provides an effective alternative approach to conventional FE methods that yields strain and stress distributions that are comparable with those generated by conventional FE methods for large deformation structural mechanics.
The primary focus of this study is to verify and to validate the stress and strain predictions generated by the IB/FE framework by comparisons with results obtained using a previously described LV model implemented within the ABAQUS FE software and by comparisons with earlier computational and experimental data sets. For the purposes of these comparisons, it is necessary to use Figure 8. Comparison of end-diastolic pressure-volume relationships obtained using the IB/FE model and the FE results of Wang et al. [15] to measurements from ex vivo human hearts [65]. Normal, healthy heart; ICM, ischemic cardiomyopathy; DCM, diopathic dilated cardiomyopathy; LVAD, hearts supported by a left ventricular assist device.
(a) (b) (c) (d) Figure 9. Comparisons of the average logarithmic fiber strain and fiber stress between the ABAQUS-based FE and IB/FE models using the modified structural stress tensor Q Q P e , with or without the pressure shift p s and incompressibility penalty term. The average stress or strain is defined as the average value over all of the selected paths.  quasi-static conditions as well as loading conditions and, where possible, material properties corresponding to those used in these earlier studies. Our results are some of the first to consider the accuracy of the stress and strain distributions generated by the IB method for realistic models of ventricular mechanics, and they indicate that the IB method can yield quantitatively accurate stresses and strains for such models. We show that improvements in the accuracy of the stress predictions can be obtained by a simple modification to the Lagrangian stress tensor that acts to reduce the magnitude of discontinuities in the Eulerian pressure field. We further show that obtaining high accuracy can require imposing the incompressibility constraint in both Lagrangian and Eulerian forms. This latter result is contrary to the prevailing approach used within the IB community, whereby it is typically assumed that it is sufficient to impose incompressibility only as an Eulerian constraint and not also as a Lagrangian constraint [21]. Although good agreement is demonstrated between the IB/FE-based model and the benchmark FE model implemented in the ABAQUS FE software, the strain and stress results generated by the current model via either approach need to be interpreted with some care because of limitations of the model: (i) the material parameters are not subject specific but rather are obtained from a porcine study [46]; (ii) the endocardial pressure loads are uniform and are no subject-specific; (iii) the heart geometry is simplified and does not include a right ventricle, the atria, the valves, and the great vessels; (iv) LV biomechanics under diastolic loading is modeled as purely passive inflation, as is commonly done in the literature [13]; and (v) experimental results from different animal models are used for strain comparison as in Figure 7.
We remark that the present model also does not account for the dynamics of the diastolic phase of the cardiac cycle. Instead, it yields stress and strain distributions that correspond to late diastole (after diastasis [66]), when the heart has been passively filled and the heart is approximately in an equilibrium state, or to ex vivo conditions in which the heart is externally pressurized, as in the experiments described in [63][64][65]. A more complete model of the diastolic phase of the cardiac cycle would include the rapid filling, slow filling, and atrial contraction phases of diastole. Such a model would require a description of cardiac dynamics. Moreover, obtaining the correct stress and strain states within the LV myocardium at the onset of diastole would seem to require a model of myocardial active tension generation and relaxation, as well as a description of both the strain state and the activation state at the onset of diastole. Although these features can all be incorporated into the present modeling framework, they are beyond the scope of this study, which focuses on verifying and validating the computational predictions of this IB/FE model under typical diastolic pressure loads.
In addition, there remain some relatively small discrepancies between the IB/FE-based and FEbased models. For example, the maximum displacement obtained from the IB/FE model is within 1% of that obtained using the FE-based model. Fiber stresses and strains along the selected transmural paths across the left ventricular free wall are also slightly different. Stress and strain are derived quantities in both the IB/FE-based and FE-based models, and it is well known that differences in displacements can lead to differences in the recovered stress or strain that are an order of magnitude larger [67]. Small differences in the computed displacements will result from different numerical treatments of the equations of structural mechanics by the IB/FE implementation and by ABAQUS and from differences in the manner in which displacement boundary conditions and pressure loading conditions are imposed in the two models. In addition, different stress recovery methods are used in the two codes. Because detailed experimental strain and stress measurements inside the LV wall are currently lacking, it is not possible to assess the actual prediction accuracy of either simulation approach, only to verify that both approaches yield comparable model predictions.
Although we only consider quasi-static conditions in this work, the IB method is well-suited for fully dynamic FSI models of the heart wall and the blood. Indeed, the IB method is primarily intended to simulate dynamics. In particular, because our IB/FE model must use a time step that is below the stability threshold associated with the explicit time stepping scheme used in our present implementation, the IB/FE model requires substantially more computational time than the comparable FE model. The IB/FE method is designed to model FSI, however, and we expect that it will be relatively straightforward to extend this model to yield subject-specific simulations of cardiac dynamics with potential clinical applications. We aim to consider the fully dynamic case in future work and to incorporate more realistic loading conditions and material models with active contraction for human ventricular tissue.

CONCLUSIONS
We have successfully applied an IB method with FE mechanics to an imaging-derived subjectspecific model of the human LV under diastolic loading conditions. The predicted stress and strain distributions agree well with our previous model, which used a conventional nonlinear FE method, as well as other experimental and computational data sets that serve to validate both the IB/FEbased and FE-based models. These results indicate that it is possible to obtain quantitatively accurate results using the IB method in conjunction with realistic experimentally constrained constitutive models and suggest that the IB method can also yield quantitative accuracy when using such constitutive models in the context of fully dynamic models of the heart.