PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of springeropenLink to Publisher's site
Cardiovascular Engineering and Technology
 
Cardiovasc Eng Technol. 2017; 8(4): 480–494.
Published online 2017 August 15. doi:  10.1007/s13239-017-0327-9
PMCID: PMC5707240

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

Abstract

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.

Keywords: Cardiac mechanics, Diastolic heart failure, Mathematical model, Myocardial strain, Pressure–volume curve, Ventricular wall stress

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. Then

x=asinhμsinνcosϕy=asinhμsinνsinϕz=acoshμcosν,
1

where (x,y,z) are Cartesian coordinates. The scale factors are gμ=gν=asinh2μ+sin2ν and gϕasinhμsinν and the Jacobian is

J=x/μx/νx/ϕy/μy/νy/ϕz/μz/νz/ϕ=acoshμsinνasinhμcosν000asinhμsinνasinhμcosν-acoshμsinν0,
2

where ϕ = 0 without loss of generality. The determinant is

|J| = a3sinhμsinν (sinh2μ + sin2ν).
3

A reference configuration is defined with a = a 0 and coordinates (μ 0, ν 0, ϕ 0), in which the LV wall occupies the region μin0 ≤ μ0 ≤ μout0 and νup ≤ ν0 ≤ π (Figs. 1a and and11b).

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 (x0y0z0) and prolate spheroidal coordinates ...

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:

|J|dμ dν dϕ = |J0|dμ0 dν0 dϕ0
4

where J 0 is the Jacobian in the reference configuration. The assumption of axisymmetry implies that dϕdϕ0. According to our previous method,27 we further assume that νν0 and obtain the mapping defined by:

a3coshμ13cosh2μ-cos2ν0-13-cos2ν0=a03coshμ013cosh2μ0-cos2ν0-a2a02coshμin013cosh2μin0-cos2ν0-a02(a0-a2)13-cos2ν0,
5

ϕϕ0a3 (cosν0 - cosνup), 
6

where aa0a1. The mapping in the μ0ν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 CFTF and the Green–Lagrange strain tensor E = ½(CI) can be computed in terms of the kinematic parameters27 (S10–15, equations labeled S refer to Supplementary Material). The dependence on cosν0 in (6) avoids a singularity in the deformation field at the apex. Equation (5) implicitly defines μ as a function of both μ0 and ν0. The dependence on ν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 μ0 (μin0 ≤ μ0 ≤ μout0), and are identified by μ0 and by their angular position ϕb (0 ≤ ϕb ≤ 2π) at the base of the heart. Position on a given fiber is parameterized by ν0 (νup ≤ ν0 ≤ π). The helical arrangement of the fibers and the variation of the helix angle through the wall are represented by setting ϕ0ϕbω(cosν0 - cosνup), where ω is a μ0-dependent wrapping parameter. As indicated in Fig. 1b, the orientation of the fibers relative to the ϕ-coordinate direction is measured either by ψ ( - π/2 < ψ ≤ π/2) or by ψ (0 ≤ ψ < π), where

ψ=π+ψ;ψ<0ψ;ψ0.
7

The angle ψ is related to ω by (S23):

cosψ=ωsinhμ0sin2ν0sinh2μ0+sin2ν0+ω2sinh2μ0sin4ν0.
8

Streeter et al.35 found that the fiber angle ψψeq at the equator of the heart (ν0π/2) varies from ψin = 85 on the inner wall to 0 at the midwall and to ψout =  - 65 on the outer wall. This is represented assuming linear variation with μ0, according to

ψeq=ψinμout-μ0μout-μin+ψoutμ0-μinμout-μin.
9

When ν0π/2, (8) gives

cosψeq=ωsinhμ0cosh2μ0+ω2sinh2μ0
10

which can be solved for the wrapping parameter,

ω=cotψeqcothμ0ifψeq0.
11

A local coordinate system (snf) is defined based on the reference fiber geometry, where f is the fiber direction, s is the transmural direction and n is perpendicular to f and s (Fig. 1b, S16–22). In the deformed geometry, the fiber paths are:

xμ0,ϕb;ν0=asinhμsinν0cos[ϕb+(ω+a3)(cosν0-cosνup)]yμ0,ϕb;ν0=asinhμsinν0sin[ϕb+(ω+a3)(cosν0-cosνup)]zμ0;ν0=acoshμcosν0
12

and the increment of arc length along the fiber is given by:

ds=a(1+(μ/ν0)2)(sinh2μ+sin2ν0)+(ω+a3)2sinh2μsin4ν0dν0,
13

where μ/ν0 is given by (S7). In the reference configuration,

ds0=a0sinh2μ0+sin2ν0+sinh2μ0sin4ν0ω2dν0.
14

The fiber stretch ratio is λfds/ds0Ls/Ls0 where L s is the current sarcomere length and L s0 is the length of the sarcomere in the reference configuration. When ψeq = 0, ω is undefined and the fibers are parameterized instead by:

xμ0,ν0;ϕ=asinhμsinν0cosϕyμ0,ν0;ϕ=asinhμsinν0sinϕzμ0,ν0=acoshμcosν0.
15

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:

σffA(t)G(Ls)(F0kavdεf/dt)(1 + kεf,ed), 
16

GLs=exp-Ls-Lsmax22Lsw2.
17

Here, F0 gives the magnitude of the stress, εf=12λf2-1 is the fiber strain, kav defines the force–velocity characteristics of the fibers and k determines the sensitivity of active force to preload according to the Frank–Starling mechanism, where εf,ed is the end-diastolic fiber strain, Lsmax is the length of the sarcomere at which maximum tension is generated and Lsw determines the width of the peak in the length–tension curve. The time-dependent fiber activation is given by

At=[sin(πt/Ta)]dwhered=1/(1+kεf,ed)when0tTa0whenTatTc,
18

where t is time after the start of contraction in a given cycle, Ta is the period of activation and Tc is the period of the cardiac cycle. The exponent d steepens the rise and fall of activation with increasing end-diastolic strain εf,ed, according to the observation that activation is prolonged with increased preload.10,34 The parameter k defines the sensitivity of the activation function to preload. Since εf,ed varies with position in the myocardium, A(t) is a function of position. The resulting Cauchy stress is σfσffef ⊗ ef where  ⊗  denotes the tensor product and ef is the basis vector for the local fiber direction. The corresponding second Piola–Kirchhoff (PK2) stress33 is Sf=λf-2σffefef.

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

Ψ = ½c1(eW - 1)
19

where

W=bffEff2+bxxEnn2+Ess2+Esn2+Ens2+bfxEfn2+Enf2+Efs2+Esf2.
20

Eij are the Green–Lagrange strain components in the fiber coordinates (snf) and c1bffbxxbfx are material parameters. The PK2 stress corresponding to σe is

Se=c1eWbxxEssbxxEsnbfxEsfbxxEnsbxxEnnbfxEnfbfxEfsbfxEfnbffEff.
21

The viscous stress σv is defined analogously to that in a viscous fluid30

σvkvm(∇v + ∇vT), 
22

where v is the material velocity field and kvm is a constant. The corresponding PK2 stress is (S32):

Sv=kvmC-1C˙C-1.
23

The components of Sf and Se are converted from fiber coordinates to prolate spheroidal coordinates using the rotation matrix (S25):

Q=1000-cosψsinψ0-sinψ-cosψ.
24

The total PK2 stress is S =  - pMC-1SfSeSv. 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,

Ω0dET:SdV0=Ωedue·TedS+Ωbdub·TbdS,
25

where Te and Tb are the boundary tractions generated by the LV pressure on the LV endocardium and the myocardial boundary at the base of the heart, dV0 is a differential volume element of the reference configuration Ω0, and Ωe and Ωb are the endocardial surface and the myocardial boundary at ν = ν up in the current configuration Ω. The incremental displacements and strains are defined by changes in the kinematic parameters, i.e., du = (u/ai)dai where u is the displacement and dE = (E/ai)dai (S33–55), and (25) gives:

Ω0ETai:SdV0=Ωeueai·TedSe+Ωbubai·TbdSbfori=1,2,3.
26

Due to axisymmetry, the LHS (internal virtual work) reduces to an integral over a rectangular region in the μ0ν0 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 ν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 ai. All terms in (26) depend nonlinearly on ai, and Sv also depends linearly on dai/dt (S68–74). Therefore, (26) leads to a system of coupled ordinary differential equations (S75) for ai(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):

Vt=Va1da1dt+Va2da2dt.
27

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 RMV and RAOV that vary depending on the pressure differences across the valves:

RMV=Rcc-Rcc-Ropmv1+e-βPA-p,
28

RAOV=Rcc-Rcc-Ropaov1+e-βp-Paov,
29

where Rcc is the resistance of a closed valve, Ropmv and Ropaov are the resistances of the open mitral and aortic valves, PA is the atrial pressure, Paov is the most proximal aortic pressure and β is a constant. The proximal aortic pressure Paov is given by

Paov=p·RAO+PAO·RAOVRAO+RAOV,
30

where RAO is the proximal aortic resistance (impedance) and PAO 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:

dPAdt=1CAPPA-PARPULM-PA-pRMV,
31

dPAOdt=1CSYSTPAOV-PAORAO-PAO-PSVRSYST,
32

Va1da1dt+Va2da2dt-PA-pRMV+p-PAOVRAOV=0,
33

where PPA is the mean pressure in the pulmonary bed, RPULM is the resistance of the pulmonary circulation and CA 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 (2833), 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 σ ff to cavity pressure p as a function of the ratio of wall volume V wall to cavity volume ratio V:

pσff=13ln1+VwallV.
34

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.

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 σ μμ, σ νν ...

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 and5b5b (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 ...

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 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 ν 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 σ 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 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 ...

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 and88 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.

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.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Animal Studies

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

Conflict of interest

All authors declare that they have no conflicts of interest.

Funding

No funding was utilized for this study.

Human Subjects

The study protocol was approved by the Institutional Review Board of the University of Nebraska Medical Center. All procedures followed were in accordance with the committee on human experimentation (institutional and national) and with the Helsinki declaration of 1975, as revised in 2000. Written informed consent was obtained from the volunteer before the study.

Footnotes

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.

References

1. Arts T, Bovendeerd PH, Prinzen FW, Reneman RS. Relation between left ventricular cavity pressure and volume and systolic fiber stress and strain in the wall. Biophys. J. 1991;59(1):93–102. doi: 10.1016/S0006-3495(91)82201-9. [PubMed] [Cross Ref]
2. Arts T, Reneman RS, Veenstra PC. A model of the mechanics of the left ventricle. Ann. Biomed. Eng. 1979;7(3–4):299–318. doi: 10.1007/BF02364118. [PubMed] [Cross Ref]
3. Augenstein KF, Cowan BR, Legrice IJ, Young AA. Estimation of cardiac hyperelastic material properties from MRI tissue tagging and diffusion tensor imaging. Med. Image Comput. Comput. Assist. Interv. 2006;9(Pt 1):628–635. [PubMed]
4. Beyar R, Sideman S. A computer study of the left ventricular performance based on fiber structure, sarcomere dynamics, and transmural electrical propagation velocity. Circ. Res. 1984;55(3):358–375. doi: 10.1161/01.RES.55.3.358. [PubMed] [Cross Ref]
5. Crisfield MA. Non-Linear Finite Element Analysis of Solids and Structures. West Sussex: Wiley; 1991.
6. Demer LL, Yin FC. Passive biaxial mechanical properties of isolated canine myocardium. J. Physiol. 1983;339:615–630. doi: 10.1113/jphysiol.1983.sp014738. [PubMed] [Cross Ref]
7. Ellwein LM, Pope SR, Xie A, Batzel JJ, Kelley CT, Olufsen MS. Patient-specific modeling of cardiovascular and respiratory dynamics during hypercapnia. Math. Biosci. 2013;241(1):56–74. doi: 10.1016/j.mbs.2012.09.003. [PMC free article] [PubMed] [Cross Ref]
8. Frenneaux M, Williams L. Ventricular-arterial and ventricular-ventricular interactions and their relevance to diastolic filling. Prog. Cardiovasc. Dis. 2007;49(4):252–262. doi: 10.1016/j.pcad.2006.08.004. [PubMed] [Cross Ref]
9. Fung YC. Biomechanics. Mechanical Properties of Living Tissues. 2. New York: Springer; 1993.
10. Gaasch WH, Carroll JD, Blaustein AS, Bing OH. Myocardial relaxation: effects of preload on the time course of isovolumetric relaxation. Circulation. 1986;73(5):1037–1041. doi: 10.1161/01.CIR.73.5.1037. [PubMed] [Cross Ref]
11. Göktepe S, Kuhl E. Electromechanics of the heart: a unified approach to the strongly coupled excitation-contraction problem. Comput. Mech. 2010;45:227–243. doi: 10.1007/s00466-009-0434-z. [Cross Ref]
12. Gorcsan J, III, Tanaka H. Echocardiographic assessment of myocardial strain. J. Am. Coll. Cardiol. 2011;58(14):1401–1413. doi: 10.1016/j.jacc.2011.06.038. [PubMed] [Cross Ref]
13. Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J. Biomech. 1995;28(10):1167–1177. doi: 10.1016/0021-9290(94)00174-3. [PubMed] [Cross Ref]
14. Guccione JM, McCulloch AD, Waldman LK. Passive material properties of intact ventricular myocardium determined from a cylindrical model. J. Biomech. Eng. 1991;113(1):42–55. doi: 10.1115/1.2894084. [PubMed] [Cross Ref]
15. Gurev V, Pathmanathan P, Fattebert JL, Wen HF, Magerlein J, Gray RA, et al. A high-resolution computational model of the deforming human heart. Biomech. Model. Mechanobiol. 2015;14(4):829–849. doi: 10.1007/s10237-014-0639-8. [PubMed] [Cross Ref]
16. Hunter PJ, McCulloch AD, ter Keurs HE. Modelling the mechanical properties of cardiac muscle. Prog. Biophys. Mol. Biol. 1998;69(2–3):289–331. doi: 10.1016/S0079-6107(98)00013-3. [PubMed] [Cross Ref]
17. Hunter PJ, Nielsen PM, Smaill BH, Legrice IJ, Hunter IW. An anatomical heart model with applications to myocardial activation and ventricular mechanics. Crit. Rev. Biomed. Eng. 1992;20(5–6):403–426. [PubMed]
18. Hunter PJ, Smaill BH. The analysis of cardiac function: a continuum approach. Prog. Biophys. Mol. Biol. 1988;52(2):101–164. doi: 10.1016/0079-6107(88)90004-1. [PubMed] [Cross Ref]
19. Janz RF, Kubert BR, Moriarty TF, Grimm AF. Deformation of the diastolic left ventricle-II. Nonlinear geometric effects. J. Biomech. 1974;7(6):509–516. doi: 10.1016/0021-9290(74)90085-2. [PubMed] [Cross Ref]
20. Kerckhoffs RC, Neal ML, Gu Q, Bassingthwaighte JB, Omens JH, McCulloch AD. Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 2007;35(1):1–18. doi: 10.1007/s10439-006-9212-7. [PMC free article] [PubMed] [Cross Ref]
21. Lafortune P, Aris R, Vazquez M, Houzeaux G. Coupled electromechanical model of the heart: parallel finite element formulation. Int. J. Numer. Method Biomed. Eng. 2012;28(1):72–86. doi: 10.1002/cnm.1494. [PubMed] [Cross Ref]
22. Landesberg A. End-systolic pressure-volume relationship and intracellular control of contraction. Am. J. Physiol. 1996;270(1 Pt 2):H338–H349. [PubMed]
23. Lumens J, Delhaas T, Kirn B, Arts T. Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction. Ann. Biomed. Eng. 2009;37(11):2234–2255. doi: 10.1007/s10439-009-9774-2. [PMC free article] [PubMed] [Cross Ref]
24. McCormick M, Nordsletten D, Kay D, Smith N. Modelling left ventricular function under assist device support. Int. J. Numer. Methods Biomed. Eng. 2011;27:1073–12095. doi: 10.1002/cnm.1428. [Cross Ref]
25. Mirsky I. Left ventricular stresses in the intact human heart. Biophys. J. 1969;9(2):189–208. doi: 10.1016/S0006-3495(69)86379-4. [PubMed] [Cross Ref]
26. Moulton MJ, Secomb TW. A low-order model for left ventricle dynamics throughout the cardiac cycle. Math. Med. Biol. 2013;30(1):45–63. doi: 10.1093/imammb/dqr024. [PubMed] [Cross Ref]
27. Moulton MJ, Secomb TW. A low-order parametric description of left ventricular kinematics. Cardiovasc. Eng. Technol. 2014;5:348–358. doi: 10.1007/s13239-014-0191-9. [Cross Ref]
28. Nordsletten DA, Niederer SA, Nash MP, Hunter PJ, Smith NP. Coupling multi-physics models to cardiac mechanics. Prog. Biophys. Mol. Biol. 2011;104(1–3):77–88. doi: 10.1016/j.pbiomolbio.2009.11.001. [PubMed] [Cross Ref]
29. Paulus WJ, Claes VA, Brutsaert DL. Physiological loading of isolated mammalian cardiac muscle. Circ. Res. 1976;39(1):42–53. doi: 10.1161/01.RES.39.1.42. [PubMed] [Cross Ref]
30. Rajagopal KR. A note on a reappraisal and generalization of the Kelvin–Voigt model. Mech. Res. Commun. 2009;36(2):232–235. doi: 10.1016/j.mechrescom.2008.09.005. [Cross Ref]
31. Regen DM. Myocardial stress equations: fiber stresses of the prolate spheroid. J. Theor. Biol. 1984;109(2):191–215. doi: 10.1016/S0022-5193(84)80003-X. [PubMed] [Cross Ref]
32. Schmid H, O’Callaghan P, Nash MP, Lin W, Legrice IJ, Smaill BH, et al. Myocardial material parameter estimation: a non-homogeneous finite element study from simple shear tests. Biomech. Model. Mechanobiol. 2008;7(3):161–173. doi: 10.1007/s10237-007-0083-0. [PubMed] [Cross Ref]
33. Spencer AJM. Continuum Mechanics. London: Longman Ltd.; 1980.
34. Spratt JA, Tyson GS, Glower DD, Davis JW, Muhlbaier LH, Olsen CO, et al. The end-systolic pressure-volume relationship in conscious dogs. Circulation. 1987;75(6):1295–1309. doi: 10.1161/01.CIR.75.6.1295. [PubMed] [Cross Ref]
35. Streeter DD, Jr, Spotnitz HM, Patel DP, Ross J, Jr, Sonnenblick EH. Fiber orientation in the canine left ventricle during diastole and systole. Circ. Res. 1969;24(3):339–347. doi: 10.1161/01.RES.24.3.339. [PubMed] [Cross Ref]
36. Suga H, Sagawa K. Mathematical interrelationship between instantaneous ventricular pressure-volume ratio and myocardial force-velocity relation. Ann. Biomed. Eng. 1972;1(2):160–181. doi: 10.1007/BF02584205. [PubMed] [Cross Ref]
37. ter Keurs HE, Rijnsburger WH, van Heuningen R, Nagelsmit MJ. Tension development and sarcomere length in rat cardiac trabeculae. Evidence of length-dependent activation. Circ. Res. 1980;46(5):703–714. doi: 10.1161/01.RES.46.5.703. [PubMed] [Cross Ref]
38. Yin FC. Ventricular wall stress. Circ. Res. 1981;49(4):829–842. doi: 10.1161/01.RES.49.4.829. [PubMed] [Cross Ref]
39. Yin FC, Chan CC, Judd RM. Compressibility of perfused passive myocardium. Am. J. Physiol. 1996;271(5 Pt 2):H1864–H1870. [PubMed]

Articles from Springer Open Choice are provided here courtesy of Springer