Determinants of biventricular cardiac function: a mathematical model study on geometry and myofiber orientation

In patient-specific mathematical models of cardiac electromechanics, usually a patient-specific geometry and a generic myofiber orientation field are used as input, upon which myocardial tissue properties are tuned to clinical data. It remains unclear to what extent deviations in myofiber orientation and geometry between model and patient influence model predictions on cardiac function. Therefore, we evaluated the sensitivity of cardiac function for geometry and myofiber orientation in a biventricular (BiV) finite element model of cardiac mechanics. Starting out from a reference geometry in which myofiber orientation had no transmural component, two new geometries were defined with either a 27 % decrease in LV short- to long-axis ratio, or a 16 % decrease of RV length, but identical LV and RV cavity and wall volumes. These variations in geometry caused differences in both local myofiber and global pump work below 6 %. Variation of fiber orientation was induced through adaptive myofiber reorientation that caused an average change in fiber orientation of ∼8∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\sim }8^\circ $$\end{document} predominantly through the formation of a component in transmural direction. Reorientation caused a considerable increase in local myofiber work (∼18%)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\sim }18\,\%)$$\end{document} and in global pump work (∼17%)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\sim }17\,\%)$$\end{document} in all three geometries, while differences between geometries were below 5 %. The findings suggest that implementing a realistic myofiber orientation is at least as important as defining a patient-specific geometry. The model for remodeling of myofiber orientation seems a useful approach to estimate myofiber orientation in the absence of accurate patient-specific information.


Introduction
In the last decade, several biventricular (BiV) finite element (FE) models have been developed as a step toward model-assisted clinical decision-making in cardiac disease (Aguado-Sierra et al. 2011;Gurev et al. 2011;Niederer et al. 2011;Sermesant et al. 2012). These models take into account the heart's geometry, myofiber orientation, and myocardial tissue properties. Whereas geometrical characteristics can be captured from imaging data and included into models through sophisticated meshing algorithms (Lamata et al. 2014), determination of local myofiber orientation and material properties is more challenging. Usually generic data from in vitro DTMRI or histological measurements (LeGrice et al. 1997;Nielsen et al. 1991) are morphed into geometrically correct models (Peyrat et al. 2007;Vadakkumpadan et al. 2012). Next, myocardial tissue properties are tuned for maximum agreement between model predicted and clinically measured cardiac function (Aguado-Sierra et al. 2011;Niederer et al. 2011 Fig. 1 The model of BiV mechanics. a Illustration of geometrical parameters, used in Table 1, in cross sections of the REF geometry. b FE mesh of simulation REF, and myofiber orientation vector in the unloaded state e f,0 quantified by helix angle α h,0 and transverse angle α t,0 using a local cardiac coordinate system { e t,0 , e l,0 , e c,0 }. c Elongated mesh (LONG) and a mesh with a higher attachment of the right to the left ventricle (RVAT). d Lumped parameter model of the circulation. The valves are modeled as ideal diodes: AV aortic valve, MV mitral valve, PV pulmonary valve, TV tricuspid valve. Both pulmonary and systemic circulation are modeled using compliances and resistances: C art , arterial compliance; C ven , venous compliance; R art , arterial resistance; R per , peripheral resistance; R ven , venous resistance; V art,0 , zero-pressure arterial volume; V ven,0 , zero-pressure venous volume. Parameter values are listed in Table 2 generic and the (unknown) patient-specific field of myofiber orientation.
The accuracy of DTMRI is on the order of ±10 • (Reese et al. 1995), even in the in vitro case (Hsu et al. 1998;Scollan et al. 1998). This accuracy seems to be sufficient for modeling cardiac electrophysiology (Vadakkumpadan et al. 2012). However, in models of left ventricular (LV) mechanics it has been shown that variations in the spatial distribution of myofiber orientations of 8 • introduced variations in myofiber stress of 10 % (Geerts-Ossevoort et al. 2003). In the same model, a drastic geometrical variation in the short-to long-axis ratio from 1:1 to 1:5 yielded similar variations in myofiber stress of 10 %. However, these results might not be valid in a BiV model, among others because of the nonsymmetrical BiV shape and the attachment of the right ventricle (RV) to the LV.
The aim of this study was to investigate the sensitivity of cardiac function to geometry and myofiber orientation in a BiV FE model of cardiac mechanics. First, local myofiber function and global pump function were compared between three geometries, which differed only by either a 27 % difference in LV short-to long-axis ratio or a 16 % difference in RV length. Next, in each of these three models the sensitivity for myofiber orientation was tested by letting the myofiber orientation change according to an adaptive myofiber reorientation model. We evaluated sensitivity to geometry and fiber orientation by assessing global function, local function, and cardiac deformation.

Geometry
In the passive stress-free state, the BiV geometry was described by two intersecting, truncated ellipsoids ( Fig. 1), defined by 8 input parameters (Table 1). Left ventricular wall volume V lvw determines the overall size, whereas dimensionless parameters f V , f sep , f R , f h , f T , f V RL , and θ A determine BiV shape. Parameter values were chosen such that the geometry for the reference simulation REF was representative for an average healthy human heart (Ho and Nihoyannopoulos 2006).

Myofiber orientation
Myofiber orientation in the unloaded state e f,0 was defined in a local cardiac coordinate system { e t,0 , e l,0 , e c,0 } by two angles (Fig. 1b). The helix angle α h,0 was defined as the angle between e c,0 and the projection of e f,0 on the circumferentiallongitudinal plane ( e c,0 , e l,0 ). The transverse angle α t,0 was defined as the angle between e c,0 and the projection of e f,0 on the circumferential-transmural plane ( e c,0 , e t,0 ). We used the bidirectional spherical linear interpolation method of Bayer et al. (2012) to ensure a smooth transition of the coordinate system and myofiber orientation in all regions, including the RV-LV attachment.
LV cavity-to-wall volume ratio in the stress-free state RV-to-LV cavity volume ratio Endocardial LV radius-to-length ratio Truncation height above equator-to-LV endocardial length ratio LV-to-RV wall thickness ratio Lowest attachment angle of RV to LV See Fig. 1a for definition of R's, Z 's, h, and θ A a The value between brackets holds for simulation LONG b The values between brackets hold for simulation RVAT During simulation of the cardiac cycle, deformation of the tissue caused a change from the reference orientation e f,0 into an actual myofiber orientation e f . In addition, e f,0 was subject to adaptation that caused a structural change in myofiber orientation in the unloaded reference state as proposed by Kroon et al. (2009): with an adaptation time constant κ that was set to four times the cardiac cycle time. e * f is the actual myofiber orientation in the deformed tissue corrected for rigid body rotation: with λ f the myofiber stretch ratio and F the deformation gradient tensor that consists of the deformation U and rigid body rotation R. In nodes on the endo-and epicardial surfaces, the myofiber vector resulting from adaptation was modified by projecting it onto the surface, to ensure that myofibers do not stick out of these surfaces.

Material properties
Myocardial tissue Cauchy stress σ was composed of a passive component σ p and an active component σ a : Constitutive relations for σ p and σ a were identical to those used by Kerckhofs et al. (2003). Passive material behavior was assumed nonlinearly elastic, transversely isotropic, and nearly incompressible, with parameter settings identical to those used by Kerckhofs et al. (2003). Active stress σ a depends on time elapsed since activation t a , sarcomere length l s , and sarcomere shortening velocity, with parameter settings identical to those used by Bovendeerd et al. (2009). Active stress development was initiated simultaneously throughout both ventricles with a cycle time of 800 ms.

Governing equations and boundary conditions
In the model, the quasi-static equations of conservation of linear momentum were solved: with ∇ the spatial gradient operator. At the base, rotation of the endocardial basal ring and movement of the whole basal plane in normal direction were prevented. The epicardial surface was assumed to be traction free while the RV and LV endocardial surfaces were uniformly subjected to a right and left ventricular pressure, p r v and p lv , respectively. During isovolumic contraction (IC) and relaxation (IR) phases of the cardiac cycle, p r v and p lv were determined such that mechanical equilibrium of the myocardial tissue was obtained at a constant end-diastolic or end-systolic RV and LV volume, respectively. During the filling and ejection phase, p r v and p lv were computed from the interaction of the RV and LV with lumped parameter models of the pulmonary and systemic circulation, respectively ( Fig. 1d) . Parameter values of the resistances and compliances were chosen such that they are representative for characteristics of the human circulation (Table 2) (Shoukas 1975).

Numerical implementation
The equilibrium equations (4) were solved numerically with a Galerkin-type FE method using 27-noded hexahedral elements with a tri-quadratic interpolation of the displacement field. The mesh consisted of 684 elements and 19803 degrees of freedom. The time step was fixed at 1 ms. It took about 2 h to compute one cardiac cycle on a single core of an Intel Core i7-2600 processor on a personal computer.

Simulations performed
The sensitivity of cardiac function to geometry was investigated using three different geometries with identical LV wall volume of 160 ml and RV wall volume of 40 ml. In the unloaded state, both LV and RV cavity volume were 60 ml. With respect to the reference geometry REF ( Fig. 1b; Table 1), parameter f R was decreased from 0.55 to 0.40 to create a more elongated geometry LONG (Fig. 1c, top). Geometry RVAT was obtained by changing θ A from 0.85π to 0.75π , thus shifting the attachment of the RV toward the base (Fig. 1c, bottom). To maintain the same RV wall volume, the LV-to-RV wall thickness ratio f T was changed from 0.33 to 0.35. In all geometries, the initial distribution of α h,0 varied nonlinearly with the transmural position from 65 • at the endocardium to −50 • at the epicardium ). The initial condition for α t,0 was set to zero. Still, a transmural component in myofiber orientation is already observed near the attachment regions (see dashed lines in septal long-axis cross sections in Fig. 4) as a result of the interpolation method (Bayer et al. 2012).
The sensitivity of cardiac function to myofiber orientation was investigated by changing myofiber orientation through adaptive reorientation. The first 12 consecutive cardiac cycles were used to approximate a hemodynamic steady state while myofiber reorientation was disabled. In the subsequent 20 cardiac cycles, myofiber orientation was adapted per node according to the adaptation model (Eq. 1). Adaptation was disabled in the final 5 cardiac cycles, during which a steady hemodynamic state was reached again.
Results from the three geometries were analyzed in terms of global pump function and local myofiber function in both the initial state (averaged over cycles 8-12) and the adapted state (averaged over cycles 33-37). Global function was quantified as cardiac output (CO), LV ejection fraction (EF), and stroke work W . Local function was quantified as natural myofiber strain during isovolumic contraction f,ic , during ejection f,ej , and during isovolumic relaxation f,ir , as well as by maximum active myofiber Cauchy stress σ f,max , and stroke work density w f calculated from the area within the active myofiber Cauchy stress-natural strain loop. Local function variables are presented as mean and standard deviation (SD) calculated from the nodes in the four LV and septal cross sections, indicated in the mesh in Fig. 4. We also evaluated myofiber angles α h,0 and α t,0 in the adapted state. Finally, we compared model computed with experimentally determined circumferential-radial shear E cr , referred to the state at begin ejection Delhaas et al. 2008;Pluijmert et al. 2014;Ubbink et al. 2006).

Results
In the initial state, only small differences in global pump function were present between geometries ( Fig. 2; Table 3). For   (Table 4): Myofiber shortening during ejection increased by 15, 9, and 11 %, σ f,max increased by 4, 4, and 5 %, and w f increased by 18, 14, and 18 %, respectively. Absolute myofiber strain during the isovolumic phases decreased. All improvements in local function were statistically significant ( p < 0.01), except for the increase of σ f,max in simulation REF( p = 0.057), LONG ( p = 0.021), and RVAT ( p = 0.013). Differences in local function were not statistically different in between geometries with adapted myofiber orientation. Figure 3 shows the mean ± SD of the differences in local function variables between corresponding nodal points in the initial and adapted state, and between corresponding nodal points in different geometries. Increase of myofiber strain during ejection (top row, middle panel) due to reorientation is significantly different from zero for all three geometries (first three values). To a lesser extent, this also holds for changes in myofiber strain during the isovolumic phases, and for w f , but not for σ f,max . In both the initial and adapted state, differences between geometries are not significantly different from zero (last six values).
In the three geometries, the change in myofiber orientation was similar not only in magnitude, but also in pattern. The helix angle α h,0 hardly changed, implying that the change in orientation took place especially in the transmural direction. Results on the transverse angle α t,0 are shown in detail in Fig. 4 for simulation REF. In the LV free wall, a gradient in α t,0 developed with negative values near the apex and positive values near the base. Largest values for α t,0 developed in the septum, especially at the RV endocardial side. In the basal region of the septum, α t,0 was positive at the posterior side, but about zero at the anterior side. To a lesser extent, asymmetry in myofiber orientation between anterior and posterior side was also seen in the LV and RV free wall. In the RV free wall, α t,0 showed similar patterns as in the LV free wall, but angles are smaller compared to the LV free wall and septum.
Temporal variations in circumferential-radial shear E cr in simulation REF are shown in Fig. 5. In the experimental results (Pluijmert et al. 2014), the total amplitude of E cr is ∼0.10 and the gradient in E cr at end ejection is apex-midbase. In the initial state in the model, the total amplitude of E cr is ∼0.35 with a base-mid-apex gradient at end ejection. In the adapted state, the total amplitude reduced to ∼0.15 and the gradient at end ejection flipped to apex-mid-base, which agrees better with the experimental data. Still, clear differences remain between patterns of E cr in model and experiment. Again, a similar trend was observed in simulations LONG and RVAT.

Discussion
We investigated the sensitivity of cardiac function for geometry and myofiber orientation in a BiV FE model of cardiac Table 4 Local function presented as (mean ± SD) natural myofiber during isovolumic contraction f,ic , during ejection f,ej , and during isovolumic relaxation f,ir , as well as maximum active myofiber Cauchy stress σ f,max , and stroke work density w f . Mean and SD were calculated from the septal and LV free wall nodes indicated in Fig. 4

Model setup
We used a parameterized description of the BiV geometry, since it enabled well-defined variations of the geometry while keeping wall and cavity volumes constant. This was important, since it is known that the ratio of cavity-to-wall volume largely determines the mean level of active wall stress (Arts et al. 1991). We also started out from a parameterized description of the myofiber orientation. While it seems tempting to vary orientation by adapting coefficients in this description, we considered such variations rather arbitrary. By using the evolution equation (1), we applied a less arbitrary, and presumably more physiologically based variation. Finally, values of the venous compliances, C s,ven and C p,ven were reduced tenfold in comparison with reported experimental data (Shoukas 1975). This choice lowered the typical time constant of the circulation model and thus decreased the num-ber of cycles needed to reach a hemodynamic steady state.
Since the same circulation model was used in all simulations, we expect that this modification did not affect the results of our sensitivity study.

Relation between myofiber orientation and cardiac function
To explain the increase in LV hemodynamic function, we consider sarcomere length in the LV and Sept nodes, indicated in Fig. 4. In model REF at end filling, sarcomere length is smaller in the adapted state (2.16 ± 0.05 µm) than in the initial state (2.18 ± 0.05 µm), in agreement with the lower end-diastolic volume in the adapted state (Fig. 2). During the isovolumic contraction (IC) phase ( f,ic ), myofiber shortening is much lower in the adapted state (0.02 ± 0.03 µm) than in the initial state (0.08±0.04 µm). Thus, sarcomere length at begin ejection is higher in the adapted state (2.14±0.07 µm) than in the initial state (2.10 ± 0.09 µm), enabling a stronger contraction because of the Frank-Starling mechanism. This causes an increased ejection pressure in the adapted state, causing an increase in arterial blood volume. Since total blood volume is constant, this leads to a decrease in venous blood volume, which causes the reduction of end-diastolic filling pressure and volume, observed in Fig. 2. The reduced myofiber shortening during IC, as caused by the change in myofiber orientation, is reflected also in a lower E cr during IC. The latter agrees better with experimental data: Although these data do not contain the IC phase, the results in Fig. 5 suggest that little E cr occurs during IC, as the signals tend to go back to zero toward the end of filling.

Remodeling of myofiber reorientation
The initial distribution for α h,0 was based on a distribution that was already optimized for homogeneity of myofiber strain during ejection in the LV (Rijcken et al. 1999). This explains why changes in α h,0 between the initial and adapted state were small. To investigate the properties of the adapted state, six additional simulations were performed, in which the adapted fiber orientation in model REF was globally changed: We increased and decreased α h,0 by 5 • , we multiplied α h,0 by 1.1 and 0.9, and we multiplied α t,0 by 1.1 and 0.9. These variations all reduced global and local function (on average 2 % reduction in both pump work and myofiber shortening during ejection) and increased standard deviations of local function parameters, reflecting a decrease in homogeneity. This indicates that the solution from the adaptation model is well defined. Experimental verification of the adapted distribution of α h,0 and α t,0 is difficult because of the lack and inaccuracy of experimental data. A difference in myofiber orientation of 8 • can hardly be detected with in vitro, let alone with in vivo DTMRI (Scollan et al. 1998;Tseng et al. 1999). The adapted distribution of the helix angle α h,0 falls within the admittedly wide range of experimental data. Model predictions of the developed α t,0 in the LV free wall and differences in α t,0 in between the anterior and posterior septum (Fig. 4) are supported by results from DTMRI studies (Geerts-Ossevoort et al. 2002;Helm et al. 2005). Since we found circumferential-radial shear strain E cr to be very sensitive to myofiber orientation, we used E cr as an indirect measure to test fiber orientation. The better match between model computed and experimentally determined E cr in the adapted state, and the finding that the adapted fiber field is well-defined, suggest that myofiber reorientation could be considered as a method to implement myofiber orientation in (geometrically correct) BiV FE models.

Toward patient-specific modeling
In recent BiV patient-specific FE models, generic data on myofiber orientation were used. By tuning model parame- ters other than myofiber orientations, good agreement was reported between model computed and experimental data on hemodynamic function (Sermesant et al. 2012). Still, as not only shown in our study but also suggested from other studies (Carapella et al. 2014), a match in hemodynamic function is no guarantee for correctly predicting deformation. Deviations in circumferential-radial and circumferentiallongitudinal shear strain from experimental observations (Carapella et al. 2014) might result from the absence of a transverse component in myofiber orientation. The model of reorientation of myofiber orientation might be used as a method to implement myofiber orientation in (geometrically correct) BiV FE models.

Study limitations
The results from our model study are dependent on material properties, in particular properties affecting shear deformation. We used a transversely isotropic constitutive model for the passive material that cannot capture the orthotropic behavior that arises from the structural organization of the myocardial tissue in sheets (LeGrice et al. 1997). It has also been shown that including active stress development perpen-dicular to the myofiber direction influences the amplitudes of the developed α t,0 , but not the gradient of α t,0 after reorientation (Pluijmert et al. 2014). The effect of sheet orientation and orthotropy needs further investigation. Results of the adapted state were analyzed after 20 cardiac cycles with adaptive myofiber reorientation. At this stage, myofiber reorientation did not yet reach a steady state: Myofiber orientation still changed by about 1 • per cardiac cycle. Because local and global function did not change anymore, we considered the situation after 20 cycles of adaptation representative for the effect of adaptation. Finally, experimental verification of our results on sensitivity of cardiac function to geometry and fiber orientation is difficult, since this would require the ability to select real hearts that differ only in either geometry or fiber orientation. Still, since we only studied a limited number of variations of geometry and fiber orientation, an extended sensitivity analysis might be performed to check our results.

Conclusion
This model study shows that a geometrical variation of 27 % in the LV short-to long-axis ratio or a variation of 16 % in RV length leads to a change in total pump work of about 5 %, whereas an average change in myofiber orientation of about 8 • causes an increase in total pump work of 11-19 % in all three geometries. Differences in total pump work between the three geometries with adapted myofiber orientation are below 5 %. These findings indicate the importance of a realistic choice of myofiber orientation. As current techniques for measurement of myofiber orientation lack sufficient accuracy, the model for remodeling of myofiber orientation seems a useful approach to set myofiber orientations in BiV FE models.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.