Simulation of Left Ventricular Dynamics Using a Low-Order Mathematical Model

The eventual goal of this study is to develop methods for estimating dynamic stresses in the left ventricle (LV) that could be used on-line in clinical settings, based on routinely available measurements. Toward this goal, a low-order theoretical model is presented, in which LV shape is represented using a small number of parameters, allowing rapid computational simulations of LV dynamics. The LV is represented as a thick-walled prolate spheroid containing helical muscle fibers with nonlinear passive and time-dependent active contractile properties. The displacement field during the cardiac cycle is described by three time-dependent parameters, using a family of volume-preserving mappings based on prolate spheroidal coordinates. Stress equilibrium is imposed in weak form and the resulting force balance equations are coupled to a lumped-parameter model of the circulation, leading to a system of differential–algebraic equations, whose numerical solution yields predictions of LV pressure and volume, together with spatial distributions of stresses and strains throughout the cardiac cycle. When static loading of the passive LV is assumed, this approach yields displacement and stress fields that closely match results from a standard finite-element approach. When dynamic motion with active contraction is simulated, substantial variations of fiber stress and strain through the myocardium are predicted. This approach allows simulations of LV dynamics that run faster than real time, and could be used to determine patient-specific parameters of LV performance on-line from clinically available measurements, with the eventual goal of real-time, patient-specific analysis of cardiac parameters. Electronic supplementary material The online version of this article (doi:10.1007/s13239-017-0327-9) contains supplementary material, which is available to authorized users.


INTRODUCTION
The pumping function of the left ventricle (LV) depends not only on the active contractile force generated in its myocardial wall, but also on the size and shape of the LV, the thickness and passive mechanical properties of the wall, and the interaction of the LV with other parts of the circulatory system. Numerous theoretical models for LV mechanics have been developed for investigating normal LV function and the pathophysiology of diseased hearts, and for eventual application in patient-specific diagnosis and treatment planning.
One general class of models is based on geometrically detailed representations of cardiac anatomy, mechanics, electrophysiology and cardiac-circulatory system interactions. 11,13,17,20,28 Typically, these models are analyzed using the finite element method, which uses a large number of degrees of freedom to represent cardiac deformation and other variables of interest. Such models have been used to estimate active and passive cardiac muscle parameters, by minimizing differences between predicted strains and experimental observations, for example from magnetic resonance imaging tagging. 3,32 The computational cost of such models is generally high, and simulation of a single cardiac cycle may take hours or longer. 15,21,24 Faster computational models for LV mechanics can be developed by restricting the number of degrees of freedom used to describe LV deformation. In this class of models, referred to here as ''low-order,'' the shape of the LV is represented approximately by a limited family of shapes, which can be specified by a small number of parameters. Generally, the shapes considered have rotational symmetry about the longitudinal axis of the LV. Examples of low-order models for LV dynamics include cylindrical models, 2,26 spherical models 23 and prolate spheroidal models. 4,25,31 A low-order model that requires a rotationally symmetric LV shape 1 is based on the assumption that muscle fiber stress and strain are approximately uniform throughout the myocardium. This model leads to the prediction that the ratio of fiber stress to LV internal pressure depends only on the ratio of LV wall volume to cavity volume, independent of LV shape. The varying elastance model, 36 in which the mechanical characteristics of the pumping LV are represented by a single scalar equation, can be considered as an extreme of the low-order approach. This model is limited in that its parameters are not directly derived from mechanical characteristics of the myocardium and it does not include an explicit representation of myocardial stress and strain.
Advances in speckle tracking echocardiography have led to the availability of data on ventricular motion and deformation with excellent temporal resolution throughout the cardiac cycle. This opens up the possibility of deducing basic muscle parameters ''on-line'' in clinical settings 12 based on ultrasound imaging. A system for patient-specific imaging and modeling, running on an echocardiography machine in a physician's office or intensive care unit, would potentially allow physicians to assess cardiac muscle properties in real time, and to predict responses to therapeutic interventions. This would involve running multiple simulations of the cardiac cycle over ranges of parameters describing muscle properties to obtain the best fit to imaging data, and would therefore require that individual simulations run much faster than real time. For this purpose, loworder models as described above are more suitable than geometrically detailed models.
To be useful for this application, a low-order model should be able to approximately represent the shape and deformation of the LV through the cardiac cycle. In this regard, a prolate spheroidal model provides a better approximation than a cylinder or a sphere, because it includes the variation in wall curvature from base to apex. According to the law of Laplace, the generation of pressure is highly dependent on wall curvature. With regard to LV deformation, three main modes of deformation can be identified, namely baseto-apex shortening, circumferential contraction and torsion, and these should all be represented in the model. The myocardium is essentially incompressible and so its motions should be restricted to volumepreserving deformations. None of the models mentioned earlier meets all these criteria. For example, previous dynamic models based on prolate spheroidal geometry 4,25 assumed a fixed relationship between the major and minor axes of the elliptical profile during deformation and therefore represent only a single mode of deformation. As a step towards developing a model for LV dynamics that meets the requirements outlined above, we previously proposed a low-order model for LV kinematics, 27 in which the LV is represented as a thickwalled, axisymmetric prolate spheroid that deforms according to a family of volume-preserving mappings defined by three time-dependent variables, and demonstrated the capability of the model to approximate the strain fields observed in the LV using speckle tracking echocardiography. A key aspect of that work was the development of volume-preserving mappings that represent three main modes of deformation, namely LV shortening, contraction and torsion, including nonlinear kinematic effects due to large strains. 19 Here, we use this representation of LV kinematics by a small number of parameters as the basis for a low-order model of LV dynamics, and show that the model can be used to predict distributions of fiber stress and through the myocardium in computational simulations that run faster than real time.

Overview
The goal of the present study is to develop a loworder model for LV dynamics using this family of mappings, based on the theory of large-deformation continuum mechanics, including nonlinear active and passive material properties, viscoelasticity, and anisotropy with respect to muscle fiber orientation. 6,9,13,18 A system of helical cardiac muscle fibers is introduced within the wall of the prolate spheroid, with fiber angle varying through the wall. The fiber angle and fiber strain are derived as functions of the kinematic parameters. The passive properties of the myocardium are represented by a Kelvin-Voigt type viscoelastic model. The elastic component is assumed to be transversely isotropic with respect to the fiber direction and is represented using a strain-energy function. The viscous component is characterized by a linear dependence on strain rate. Active fiber stress is dependent on a prescribed time-dependent activation function, and takes into account the length-tension and force-velocity characteristics of cardiac muscle, and the dependence of force generation on end-diastolic strain (Frank-Starling effect). The equations of stress equilibrium are solved approximately in weak form. The resulting equations are coupled to a lumped-parameter model of the circulatory system. This allows the dynamic behavior of the LV and its interaction with the circulatory system to be represented by a system of differential-algebraic equations, which can be solved rapidly by standard methods.

Coordinate Systems
Prolate spheroidal coordinates (l,m,/) are used with a time-dependent interfocal distance 2a. Then where (x,y,z) are Cartesian coordinates. The scale factors are g l ¼ g m ¼ a ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi sinh 2 l þ sin 2 m q and g / ¼ a sinh l sin m and the Jacobian is J ¼ @x=@l @x=@m @x=@/ @y=@l @y=@m @y=@/ @z=@l @z=@m @z=@/ where a ¼ a 0 þ a 1 . The mapping in the l 0 À m 0 plane is illustrated in Fig. 1c. The first kinematic parameter gives lengthening (a 1 > 0) or shortening (a 1 < 0) of the LV at almost constant volume. The second parameter gives contraction (a 2 > 0) or expansion (a 2 < 0) (Fig. 1d). The third parameter gives torsional deformation that increases from base to apex and is clockwise (a 3 < 0) or counterclockwise (a 3 > 0) as viewed from the apex. The deformation gradient tensor F, the right Cauchy-Green strain tensor C ¼ F T F and the Green-Lagrange strain tensor E ¼ 1 2 C À I ð Þ can be computed in terms of the kinematic parameters 27 (S10-15, equations labeled S refer to Supplementary Material). The dependence on cos m 0 in (6) avoids a singularity in the deformation field at the apex. Equation (5) implicitly defines l as a function of both l 0 and m 0 . The dependence on m 0 implies that the deformed inner and outer LV surfaces are not generally spheroidal.

Muscle Fiber Geometry
In the reference configuration, muscle fibers are assumed to lie in surfaces of constant l 0 ðl in0 l 0 l out0 Þ, and are identified by l 0 and by their angular position / b ð0 / b 2pÞ at the base of the heart. Position on a given fiber is parameterized by m 0 ðm up m 0 pÞ. The helical arrangement of the fibers and the variation of the helix angle through the wall are represented by setting / 0 ¼ / b þ xðcos m 0 À cos m up Þ, where x is a l 0 -dependent wrapping parameter. As indicated in Fig. 1b, the orientation of the fibers relative to the /-coordinate direction is measured either by w ðÀp=2<w p=2Þ or by w Ã ð0 w Ã <pÞ, where The angle w Ã is related to x by (S23): Streeter et al. 35 found that the fiber angle w ¼ w eq at the equator of the heart ðm 0 ¼ p=2Þ varies from w in ¼ 85 on the inner wall to 0 at the midwall and to w out ¼ À65 on the outer wall. This is represented assuming linear variation with l 0 , according to When m 0 ¼ p=2, (8) gives which can be solved for the wrapping parameter, A local coordinate system ðs; n; fÞ is defined based on the reference fiber geometry, where f is the fiber direction, s is the transmural direction and n is per-pendicular to f and s (Fig. 1b, S16-22). In the deformed geometry, the fiber paths are: x 0 FIGURE 1. Coordinates and variables used to describe the geometry and deformation of the left ventricle (LV). (a) Section of reference LV shape in plane containing axis of rotational symmetry, showing Cartesian coordinates ðx 0 ; y 0 ; z 0 Þ and prolate spheroidal coordinates ðl 0 ; m 0 ; / 0 Þ. The initial wall shape is characterized by interfocal distance ð2a 0 Þ, inner boundary l 0 ¼ l in0 , outer boundary l 0 ¼ l out0 and upper edge m 0 ¼ m up . (b) Three-dimensional representation of reference configuration, showing Cartesian coordinates ðx 0 ; y 0 ; z 0 Þ and prolate spheroidal coordinates ðl 0 ; m 0 ; / 0 Þ. The local fiber coordinates ðs; n; f Þ and the fiber angles w and w Ã are shown for two helical fibers, one epicardial and one endocardial. For the epicardial fiber, w < 0, and w Ã ¼ w þ p. For the endocardial fiber, w > 0 and w Ã ¼ w. (c) Schematic diagram of mapping from reference to deformed LV shape. Points in the l 0 À m 0 plane are mapped to new l values with no change in m. In the physical x À z plane, the LV shape is deformed. (d) Effects of variations in parameters a 1 and a 2 on LV shape. In each case, the reference shape is shown by dashed curves and the deformed shape is indicated by the shaded region. and the increment of arc length along the fiber is given by: where @l=@m 0 is given by (S7). In the reference configuration, The fiber stretch ratio is where L s is the current sarcomere length and L s0 is the length of the sarcomere in the reference configuration. When w eq ¼ 0, x is undefined and the fibers are parameterized instead by:

Active Fiber Stress
The Cauchy stress generated by active fiber contraction is uniaxial in the fiber direction and experimentally has a well-defined zero tension sarcomere length and a length at which maximum tension is generated. 37 We approximate it by: Here, F 0 gives the magnitude of the stress, is the fiber strain, k av defines the forcevelocity characteristics of the fibers and k determines the sensitivity of active force to preload according to the Frank-Starling mechanism, where e f;ed is the enddiastolic fiber strain, L smax is the length of the sarcomere at which maximum tension is generated and L sw determines the width of the peak in the lengthtension curve. The time-dependent fiber activation is given by where t is time after the start of contraction in a given cycle, T a is the period of activation and T c is the period of the cardiac cycle. The exponent d steepens the rise and fall of activation with increasing end-diastolic strain e f;ed , according to the observation that activation is prolonged with increased preload. 10,34 The parameter k 0 defines the sensitivity of the activation function to preload. Since e f;ed varies with position in the myocardium, A t ð Þ is a function of position. The resulting Cauchy stress is r f ¼ r ff e f e f where denotes the tensor product and e f is the basis vector for the local fiber direction. The corresponding second Piola-Kirchhoff (PK2) stress 33

Passive Material Properties of Myocardium
The passive myocardium is represented as an incompressible anisotropic viscoelastic solid, 6,39 using a Kelvin solid model for viscoelasticity. The Cauchy stress is Àp M I þ r f þ r e þ r v , where Àp M I is the reaction stress due to the constraint of incompressibility 30 and the corresponding PK2 stress is Àp M C À1 . 33 The elastic stress r e is assumed to be transversely isotropic with respect to the fibers 14 and given by a strain-energy function 9,13 E ij are the Green-Lagrange strain components in the fiber coordinates ðs; n; fÞ and c 1 ; b ff ; b xx ; b fx are material parameters. The PK2 stress corresponding to r e is The viscous stress r v is defined analogously to that in a viscous fluid 30 where v is the material velocity field and k vm is a constant. The corresponding PK2 stress is (S32): The components of S f and S e are converted from fiber coordinates to prolate spheroidal coordinates using the rotation matrix (S25): The total PK2 stress is Because the mappings are volume-conserving, the reaction stress does no work and can be omitted from the equilibrium equations.

Equilibrium Equations
The equations of force equilibrium are expressed in weak form. 5 The total virtual work done by an incremental displacement is zero,, ZZZ where T e and T b are the boundary tractions generated by the LV pressure on the LV endocardium and the myocardial boundary at the base of the heart, dV 0 is a differential volume element of the reference configuration X 0 , and @X e and @X b are the endocardial surface and the myocardial boundary at m = m up in the current configuration X. The incremental displacements and strains are defined by changes in the kinematic parameters, i.e., du ¼ ð@u=@a i Þda i where u is the displacement and dE ¼ ð@E=@a i Þda i (S33-55), and (25) gives: ZZZ Due to axisymmetry, the LHS (internal virtual work) reduces to an integral over a rectangular region in the l 0 À m 0 plane (Fig. 1c), which is evaluated using 2D Simpson's rule on a regularly spaced 9 9 11 rectangular mesh. The first term of the RHS (external virtual work) reduces to a line integral over m 0 (S65). The second term can be computed as the virtual work done on the upper boundary of the cavity (S67). These quantities are calculated at each time point using the current values of a i . All terms in (26) depend nonlinearly on a i , and S v also depends linearly on da i =dt (S68-74). Therefore, (26) leads to a system of coupled ordinary differential equations (S75) for a i ðtÞ.

Ventricular Volume
The LV cavity volume VðtÞ is computed as a solid of revolution about the z axis and its time derivative is (S76-78):

Model of Circulatory System
The LV is coupled to a lumped-parameter model of the circulation (Fig. 2), including pulmonary vascular bed, left atrium, aorta and systemic vascular bed. Effects of LV inertia have been shown to be negligible and are not included. 7,29 The mitral and aortic valves are represented by resistances R MV and R AOV that vary depending on the pressure differences across the valves: where R cc is the resistance of a closed valve, R op m v and R op a ov are the resistances of the open mitral and aortic valves, P A is the atrial pressure, P aov is the most proximal aortic pressure and b is a constant. The proximal aortic pressure P aov is given by where R AO is the proximal aortic resistance (impedance) and P AO is the aortic pressure after the impedance. Then the differential equations describing the interaction of the LV mechanics with the pulmonary, atrial and systemic circulation are given by: @V @a 1 where P PA is the mean pressure in the pulmonary bed, R PULM is the resistance of the pulmonary circulation and C A is the compliance of the atrium.

Numerical Solution Method
The force balance Eq. (26), together with (28-33), make up a system of differential-algebraic equations. The solution is obtained through a custom algorithm that solves the time integration using second order Runge-Kutta and solves for the solution vector at each time point using the Newton-Raphson method (S76-80). The resulting system is stiff during rapid valve resistance changes (valve openings and closings) and the time step was varied accordingly.

Echocardiography Image Acquisition and Establishment of the Reference Configuration
A normal volunteer was studied with echocardiography. A standard echocardiography protocol was employed as described by Moulton and Secomb. 27 Three long-axis (AP4, AP3 and AP2 views) echocardiography images obtained from an end-systolic echo frame were used to establish a reference configuration of the LV. A prolate spheroidal reference shape was obtained by minimizing the sum of squared deviations between the fitted points on the endocardial and epi-cardial surfaces and corresponding observed points in three long-axis images. 27 The echo frame with lowest LV volume was chosen as the reference frame, given that an unloaded reference configuration is not identified by clinical imaging data.

Construction of Models to Test Effects of LV Shape
Arts et al. 1 derived an equation for the ratio of fiber stress r ff to cavity pressure p as a function of the ratio of wall volume V wall to cavity volume ratio V: This derivation applies to any rotationally symmetric shape, under the assumption that the stress at each point in the wall is the sum of a hydrostatic pressure and a uniform unidirectional fiber stress. To test the predictions of this equation, we constructed three different LV models with identical wall volumes and cavity volumes, as follows: normal, with the shape derived from an image of a normal volunteer in early diastole; spherical, with a more spherical shape; ellipsoidal, with a more prolate shape. The three shapes are shown in Fig. 3 and model parameters are given in Table 1.

Comparison with Finite Element Solution
Finite element (FE) simulations were performed using FlexPDE Version 3.11 (PDE Solutions Inc., Spokane Valley, WA 99206) for passive prolate spheroidal shells with nonlinear transversely isotropic material properties, statically loaded by internal pressure, and compared to corresponding simulations with the same loads, material properties and reference shape using the low-order method. A reference configuration 2 cm Normal Ellipsoidal Spherical FIGURE 3. LV wall shapes used to examine effects of shape on stress distribution. All three shapes have the same cavity volume, 59.8714 cm 3 , and wall volume, 158.112 cm 3 . Parameter values are given in Table 1.

Model Pressures, Volumes and Strains for Several Cardiac Cycles
The model was solved over the period of twelve cardiac cycles, for the model with normal shape and parameters given in Table 1. Solution time was 7.1 s, i.e., 0.6 s per cardiac cycle. End diastolic and endsystolic LV geometry are shown in Figs. 5a and 5b (for the twelfth cardiac cycle) along with reference fibers (solid lines) on the endocardium and epicardium  superimposed on the deformed configurations (dashed lines). Corresponding pressure-volume curves are shown in Fig. 5c. Figure 5d shows the time course of activation, LV and aortic pressures, LV volume, and strain components over the first 6 cycles of the simulation. Periodic behavior is achieved after 4-5 cycles.

Effects of LV Ellipticity on Pressure Generation and Wall Stresses
For each of the LV models with different shapes (spherical, normal, ellipsoidal) shown in Fig. 3, model solutions were obtained with the same passive and active material properties, viscous parameters and circulatory model parameters. Figure 6 shows the pressure-volume loops for the three shapes, which agree very closely, despite the different LV shapes. The stroke work values for the three models agree within about 1% (Table 1). Diastolic filling commences slightly before the end of muscle activation, resulting in a brief period during which pressure declines while volume is increasing.
By integrating over time, (25) can be used to calculate the work done during computed motions of the LV. For the normal geometry, the internal work done per cycle (averaged over twelve cardiac cycles) is 0.915 J/cycle, representing the sum of the fiber work (1.004 J/cycle) and the viscous dissipation in the myocardium (À0.089 J/cycle). The net change in elastic energy is zero over each cycle. The stroke work, computed as R pdV, is 0.900 J/cycle. Corresponding results for the spherical and ellipsoidal geometries are given in Table 1. In theory, the net internal work should equal the stroke work. The differences found here (respectively 1.6, 0.4 and 4.7%) are attributed to the approximate estimate (S67) of the work done at the base of the heart. The analysis was repeated for three corresponding models with the same internal and wall volumes, but with m up = p/2, so that there is no z displacement at the base. In that case, stroke work and net internal work agreed within 0.001%. The deviatoric components of the Cauchy stresses r // , r mm , r /m and r ff for the three shapes at peak activation are shown in Fig. 7. The spatial patterns and magnitudes of the peak systolic stresses are similar for all three LV shapes. The strong variations of circumferential and longitudinal stresses through the wall reflect the variations in fiber orientation. Fiber stress varies across the wall, with a maximum near the midwall, where fibers are nearly circumferential, near the LV equator. Towards the apex, fiber stresses are lower and the maximum occurs at the endocardium. Figure 8 shows the predicted time-dependent variation of deviatoric stress in the fiber direction at three locations through the wall (endocardium, midwall and epicardium) at m 0 = 3p/4, for the three LV shapes. The fiber stresses predicted by the model of Arts et al. 1 are also shown, as computed using (34) with the cavity pressure from the current model.

Sensitivity of Results to Reference Parameter Values
and Cavity-to-Wall Volume Ratio The predictions of these simulations were tested by varying passive muscle, active muscle and circulatory reference parameter values by ±20%. Close agreement of pressure-volume curves and stroke work and qualitative agreement of fiber stress distributions for the three geometries were found throughout the parameter range. The relationship between the Arts et al. 1 model and the current results was found to be similar throughout the range of parameters examined, and also when models with cavity to wall volume ratios of 0.2 and 0.6 were examined.

DISCUSSION
The overall goal of this work is to develop a fast computational model that links the kinematics of the beating LV with the active and passive properties of the myocardial tissue. To address this goal, we developed a model with a small number of degrees of freedom. This approach can be considered as intermediate in level of detail between the simplicity of lumped approaches such as the varying elastance model and the complexity of finite-element approaches with high spatial resolution. As such, it combines an explicit description of LV geometry and deformation with low computational cost.
The present model has several distinctive features. (i) The passive properties of the myocardium are represented using a Kelvin-Voigt viscoelastic model, with elastic properties derived from a transversely isotropic strain-energy function previously used to model the heart. 13 (ii) Nonlinear (finite deformation) kinematics and material properties are utilized. (iii) Active contraction is represented by a model of muscle that incorporates length-tension and force-velocity properties and where force generation depends on end-diastolic strain according to the Frank-Starling law. Muscle activation is a prescribed function of time, such as could be derived from an electrophysiological model. (iv) The LV is coupled to a lumped-parameter model of the circulatory system that includes opening and closing of cardiac valves. (v) The model predicts LV strains, volume, stresses and pressure throughout the cardiac cycle. (vi) Because it is a low-order model, it is executed faster than real time, with potential applications for on-line assessment of muscle properties using clinically available data including echocardiography.
The low-order model represents a compromise between physiological realism and computational speed. Only the LV is modeled, and effects of interactions with the right ventricle (RV), 8 either by direct mechanical interaction or through interactions with the pulmonary circulation, are neglected. The LV is assumed to be axisymmetric, and effects that depend on angular position, such as LV-RV interactions or localized myocardial damage, are not included. Such effects could potentially be incorporated by adding additional degrees of freedom to the mapping functions. A phenomenological model for muscle contraction is used. Incorporation of a model for the biophysics of sarcomere function 16,22 could yield more insight into the interplay between fiber contraction and relaxation and ventricular performance.
In the case of passive expansion of a pressurized spheroid, close agreement was found between the loworder model and finite element simulations with many degrees of freedom (Fig. 4). This finding suggests that a low-order approach can faithfully represent important aspects of LV mechanics.
The results in Fig. 5 illustrate the ability of the model to predict a range of parameters and behaviors that are of primary interest with regard to LV dynamics, including LV pressure and internal volume, components of stress and strain, P-V loops, and three-dimensional LV motions including contraction, shortening and torsion. This computation ran faster than real time and has potential to run much faster.
An important contribution to the understanding of the relationship between myocardial stress and LV pressure was made by Arts et al. 1 Assuming that fiber stress is uniform throughout the wall, they derived Eq. (34) giving the ratio of LV pressure to fiber stress to as a function of the ratio of wall volume to cavity volume, for any axisymmetric shape and independent of fiber orientation. This equation implies that LV shape does not have a large effect on pressure gener-ation for given fiber stress and given cavity and wall volumes. Our results are consistent with this implication. The pressure-volume loops for three LV models with matched volumes but different shapes (Fig. 3) agree closely (Fig. 6), with only 1% variation in stroke work between models (Table 1).
In contrast to the model of Arts et al., 1 the present model allows estimation of variations in fiber stress and strain through the LV wall. Figure 7 shows substantial spatial variations in fiber stress at peak systolic contraction. Figure 8 shows that the temporal variation of fiber stress, although similar in form to that predicted by the model of Arts et al., varies significantly through the wall. In the model of Arts et al., fiber stress is assumed to be uniform from endocardium to epicardium. Their model does not include passive stiffness, and does not allow prediction of diastolic stresses. Also, it does not explicitly describe torsional stresses or deformation. The results in Figs. 7 and 8 are based on an assumed reference configuration that was taken to be close to the endsystolic echo geometry. Initially, all sarcomere lengths are equal at the reference shape with value of 1.82 lm. As the heart expands and contracts, the fiber strain is spatially variable, resulting in a non-uniform distribution of fiber stress. Even if a different reference shape was assumed, there would still be transmural variations in fiber strain and therefore variations of stress, in contrast to the Arts et al. model (Table 2). Ventricular wall stress is an important determinant of muscle metabolism and growth and remodeling in normal and pathologic conditions, 38 but it is not measurable in the beating heart. Key aspects of muscle mechanics such as force-velocity and forcelength relationships are measurable in isolated muscle specimens but not in the intact LV. The model developed here provides a basis for noninvasive assessment of wall stress and active force development in the LV. Predictions of the pattern and magnitude of longitudinal, circumferential, torsional and fiber stresses, as shown in Fig. 7, have potential applications for investigating the mechanisms of normal and pathological remodeling of the myocardium, and for estimating the spatial distribution (transmural and base to apex) of metabolic needs in the LV.
A goal of this work is to utilize observed LV deformations from speckle tracking echocardiography, in conjunction with the dynamic model, to identify parameters that define active and passive muscle properties. Such parameter identification will require repeated simulations of LV motion and deformation using the model, with iterative adjustment of the parameters for best fit between predicted  and observed motions. The rapid computational speed of the present model makes this approach feasible. Alternatively, low-order model parameter estimation could serve as a preliminary step for finiteelement model parameter estimation, narrowing the parameter space and thereby reducing computational needs. The model presented here assumes axisymmetric LV geometry, but the low-order approach can be extended to the analysis of non-axisymmetric geometries by including additional modes of deformation.

ELECTRONIC SUPPLEMENTARY MATERIAL
The online version of this article (doi: 10.1007/s13239-017-0327-9) contains supplementary material, which is available to authorized users.

ANIMAL STUDIES
This article does not contain any studies with animals performed by any of the authors.