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.

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 models23 and prolate spheroidal models.4,25,31 A low-order model that requires a rotationally symmetric LV shape1 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 settings12 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, low-order 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 base-to-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 volume-preserving deformations. None of the models mentioned earlier meets all these criteria. For example, previous dynamic models based on prolate spheroidal geometry4,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 thick-walled, 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.

Methods

Overview

The goal of the present study is to develop a low-order 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 (μ,ν,ϕ) are used with a time-dependent interfocal distance 2a. Thenwhere (x,y,z) are Cartesian coordinates. The scale factors are and and the Jacobian iswhere ϕ = 0 without loss of generality. The determinant is

A reference configuration is defined with a = a 0 and coordinates (μ 0, ν 0, ϕ 0), in which the LV wall occupies the region and (Figs. 1a and 1b).

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 and prolate spheroidal coordinates . The initial wall shape is characterized by interfocal distance , inner boundary , outer boundary and upper edge . (b) Three-dimensional representation of reference configuration, showing Cartesian coordinates and prolate spheroidal coordinates . The local fiber coordinates and the fiber angles and are shown for two helical fibers, one epicardial and one endocardial. For the epicardial fiber,  < 0, and . For the endocardial fiber,  > 0 and . (c) Schematic diagram of mapping from reference to deformed LV shape. Points in the plane are mapped to new values with no change in . In the physical 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.

Volume-Preserving Deformation

The deformation of the LV wall is represented using a family of mappings (a 0;μ 0,ν 0,ϕ 0) → (a;μ,ν,ϕ) that approximates the main modes of deformation. Volume preservation is ensured if volume elements in the initial and final configurations are equal:where J 0 is the Jacobian in the reference configuration. The assumption of axisymmetry implies that . According to our previous method,27 we further assume that and obtain the mapping defined by: where . The mapping in the 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 and the Green–Lagrange strain tensor can be computed in terms of the kinematic parameters27 (S10–15, equations labeled S refer to Supplementary Material). The dependence on in (6) avoids a singularity in the deformation field at the apex. Equation (5) implicitly defines as a function of both and . The dependence on 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 , and are identified by and by their angular position at the base of the heart. Position on a given fiber is parameterized by . The helical arrangement of the fibers and the variation of the helix angle through the wall are represented by setting , where is a -dependent wrapping parameter. As indicated in Fig. 1b, the orientation of the fibers relative to the -coordinate direction is measured either by or by , where

The angle is related to by (S23):

Streeter et al.35 found that the fiber angle at the equator of the heart varies from on the inner wall to at the midwall and to on the outer wall. This is represented assuming linear variation with , according to

When , (8) giveswhich can be solved for the wrapping parameter,

A local coordinate system is defined based on the reference fiber geometry, where is the fiber direction, is the transmural direction and is perpendicular to and (Fig. 1b, S16–22). In the deformed geometry, the fiber paths are:and the increment of arc length along the fiber is given by:where is given by (S7). In the reference configuration,

The fiber stretch ratio is where L is the current sarcomere length and L is the length of the sarcomere in the reference configuration. When , 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, gives the magnitude of the stress, is the fiber strain, defines the force–velocity characteristics of the fibers and determines the sensitivity of active force to preload according to the Frank–Starling mechanism, where is the end-diastolic fiber strain, is the length of the sarcomere at which maximum tension is generated and determines the width of the peak in the length–tension curve. The time-dependent fiber activation is given bywhere is time after the start of contraction in a given cycle, is the period of activation and is the period of the cardiac cycle. The exponent steepens the rise and fall of activation with increasing end-diastolic strain , according to the observation that activation is prolonged with increased preload.10,34 The parameter defines the sensitivity of the activation function to preload. Since varies with position in the myocardium, is a function of position. The resulting Cauchy stress is where denotes the tensor product and is the basis vector for the local fiber direction. The corresponding second Piola–Kirchhoff (PK2) stress33 is .

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 , where is the reaction stress due to the constraint of incompressibility30 and the corresponding PK2 stress is .33 The elastic stress is assumed to be transversely isotropic with respect to the fibers14 and given by a strain-energy function9,13 where

are the Green–Lagrange strain components in the fiber coordinates and are material parameters. The PK2 stress corresponding to is

The viscous stress is defined analogously to that in a viscous fluid30 where v is the material velocity field and is a constant. The corresponding PK2 stress is (S32):

The components of and 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,where and are the boundary tractions generated by the LV pressure on the LV endocardium and the myocardial boundary at the base of the heart, is a differential volume element of the reference configuration , and and are the endocardial surface and the myocardial boundary at ν = ν in the current configuration . The incremental displacements and strains are defined by changes in the kinematic parameters, i.e., where u is the displacement and (S33–55), and (25) gives:

Due to axisymmetry, the LHS (internal virtual work) reduces to an integral over a rectangular region in the plane (Fig. 1c), which is evaluated using 2D Simpson’s rule on a regularly spaced 9 × 11 rectangular mesh. The first term of the RHS (external virtual work) reduces to a line integral over (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 . All terms in (26) depend nonlinearly on , and also depends linearly on (S68–74). Therefore, (26) leads to a system of coupled ordinary differential equations (S75) for .

Ventricular Volume

The LV cavity volume is computed as a solid of revolution about the 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 and that vary depending on the pressure differences across the valves: where is the resistance of a closed valve, and are the resistances of the open mitral and aortic valves, is the atrial pressure, is the most proximal aortic pressure and is a constant. The proximal aortic pressure is given bywhere is the proximal aortic resistance (impedance) and 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: where is the mean pressure in the pulmonary bed, is the resistance of the pulmonary circulation and is the compliance of the atrium.

Figure 2

Schematic diagram of circulatory model indicating resistances, pressures and compliances.

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 epicardial 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 σ to cavity pressure p as a function of the ratio of wall volume V 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.

Figure 3

LV wall shapes used to examine effects of shape on stress distribution. All three shapes have the same cavity volume, 59.8714 cm3, and wall volume, 158.112 cm3. Parameter values are given in Table 1.

Table 1

Geometric parameters and work done for the three LV shapes. Work per cardiac cycle is averaged over twelve cycles.

	Spherical	Normal	Ellipsoidal
a_0 (cm)	2.7307	4.6307	6.7307
μ_in	0.7470308	0.39521852	0.235
μ_out	1.1173771	0.6719001	0.426
ν_up	1.1112	1.1112	1.1112
V (LV volume) (cm3)	59.8714	59.8714	59.8714
V_Wall (LV wall volume) (cm3)	158.112	158.112	158.112
Fiber work/cycle (J)	1.039	1.004	0.982
Viscous work/cycle (J)	−0.105	−0.089	−0.082
Internal work/cycle (J)	0.933	0.915	0.900
External work/cycle (J)	0.933	0.915	0.900
Stroke work/cycle (J)	0.890	0.900	0.897

Results

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 was used with initial parameters a 0 = 5.0, μ in = 0.45, μ out = 0.62, ν up = π/2, yielding a cavity volume of 64.4 cm3, and material properties c 1 = 2, b ff = 4, b xx = 2, b fx = 8 in (19–20), and loaded with internal pressure 1 kPa. In the axisymmetric FE simulations, the two-dimensional domain initially contained 38 triangular elements. The mesh was refined to 183 and then 375 elements to test convergence. The current model and FE solutions are compared in Fig. 4. Displacements of the low-order model agree closely with the finite element solution with small deviations near the apex. Distributions of deviatoric stress components generally agree closely, with differences appearing most prominently in areas of high curvature. The nonlinear FE model required 5.6 h to solve and the current model required 4 s to reach an equilibrium solution on a Dell Latitude E6420 Intel Core i5-2520 M with 2.5 Ghz processor with 4 GB RAM (Dell Computer Corporation, Austin TX).

Figure 4

Comparison of current model and finite element solutions. (a) Undeformed configuration and deformed configurations computed by current model and by FE method. (b) Deviatoric normal stress components σ , σ and σ computed by current model and by FE method, and magnitude of difference. Color bar scales in kPa.

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

Figure 5

Examples of model results. (a) End diastolic configuration. (b) End systolic configuration. Reference endocardial and epicardial fiber curves are shown. LV torsion and twist of the fibers are demonstrated (dashed lines). (c) Pressure–volume curves for twelve cycles. (d). Time course of assumed activation, LV and aortic pressures, LV volume, and strain components over 6 cycles. The strain components E (circumferential), E (longitudinal) and E (torsional) are computed at a point in the mid-wall where ν 0 = 3π/4.

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.

Figure 6

Pressure–volume loops for three different reference shapes with same wall and cavity volumes. Solid line: normal. Short dashed line: spherical. Long dashed line: ellipsoidal.

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 , 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 ν up = π/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 σ , σ , σ and σ 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 mid-wall, 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 ν 0 = 3π/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.

Figure 7

Deviatoric wall stress contours at peak systole for three different reference shapes with same wall and cavity volumes. For each stress component, shapes are (left to right) spherical (S), normal (N) and ellipsoidal (E). Color bar scales in kPa.

Figure 8

Variation of deviatoric normal stress in the fiber direction through the cardiac cycle, computed at points at the endocardium, mid-wall and epicardium where ν 0 = 3π/4. Heavy curves show predictions of model by Arts et al.,1 based on uniform fiber stress throughout myocardium and the computed cavity pressures.

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 function16,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 low-order 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 generation 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 end-systolic echo geometry. Initially, all sarcomere lengths are equal at the reference shape with value of 1.82 μm. 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).

Table 2

Parameters and variables of the model.

Parameter	Value	Description
a,a_1,a_2,a_3		Kinematic parameters defining deformation field
b_ff	5	Passive elastic material coefficient
b_fx	7	Passive elastic material coefficient
b_xx	3	Passive elastic material coefficient
c_1	1	Passive nonlinear material coefficient (kPa)
C_A	20	Compliance of atrium (cm3/kPa)
C_SYST	10	Compliance of aorta (cm3/kPa)
F		Deformation gradient tensor
F_0	300	Active force generation parameter (kPa)
g_μ,g_ν,g_ϕ		Scale factors for ellipsoidal coordinates
J		Jacobian of mapping from prolate to Cartesian coordinates
k	0	Sensitivity of active force to end-diastolic fiber strain
k^′	1	Controls shape of activation function
k_av	0.5	Constant controlling amount of viscous fiber stress
k_vm	0.025	Constant reflecting the amount of viscous matrix stress
L_s0	1.82	Length of sarcomere at zero pressure (μm)
L_s		Dynamic sarcomere length (μm)
L_smax	2.23	Length of sarcomere (μm) where maximum tension is generated
L_sw	0.2	Width of exponential length–tension relationship (μm)
p		LV pressure (kPa), variable
P_A		Atrial pressure (kPa), variable
P_AO		Pressure in the proximal aorta (kPa), variable
P_AOV		Aortic pressure, just beyond aortic valve (kPa), variable
P_SV	6.0	Pressure in the systemic veins (kPa)
P_PA	3.0	Pressure in the pulmonary circulation (kPa)
R_PULM	0.023	Resistance in pulmonary circulation (kPa s/cm3)
R_MV		Resistance of mitral valve (open) (kPa s/cm3), variable
R_OP_MV	0.0008	Resistance of fully open mitral valve (kPa s/cm3)
R_AOV		Resistance of aortic valve (open) (kPa s/cm3), variable
R_OP_AOV	0.0005	Resistance of fully open aortic valve (kPa s/cm3)
R_AO	0.001	Proximal aortic resistance, constant (kPa s/cm3)
R_SYST	0.11	Resistance of systemic circulation (kPa s/cm3)
S,S_f,S_e,S_v		PK2 stress tensors: total, fiber, elastic, viscous
T_b,T_e		Boundary tractions: base of LV, endocardium
Ω_0,Ω		Domain of the LV myocardium: reference, deformed configuration
E		Green strain tensor
u		Boundary displacements for virtual work equations
T_a	0.4	Period of ventricular activation (s)
T_c	1	Period of cardiac cycle (s)
μ,μ_0		Coordinate for position through wall: deformed, reference
ν_0		Base-to-apex coordinate
ϕ,ϕ_0		Coordinate for circumferential position: deformed, reference
ψ		Local fiber angle
ψ^∗		Quadrant adjusted local fiber angle
ω		Fiber wrapping parameter
ds_0,ds		Fiber arc length: reference, deformed configuration

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 force–length 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 finite-element 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.

Below is the link to the electronic supplementary material.

Supplementary material 1 (DOCX 287 kb)