help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Supplemental Material
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Watanabe, H.
Right arrow Articles by Hisada, T.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Watanabe, H.
Right arrow Articles by Hisada, T.
Biophysical Journal 87:2074-2085 (2004)
© 2004 The Biophysical Society

Multiphysics Simulation of Left Ventricular Filling Dynamics Using Fluid-Structure Interaction Finite Element Method

Hiroshi Watanabe, Seiryo Sugiura, Hidenobu Kafuku and Toshiaki Hisada

Computational Biomechanics Division, Institute of Environmental Studies, Graduate School of Frontier Sciences, The University of Tokyo, Bunkyo-ku Tokyo 113-0033, Japan

Correspondence: Address reprint requests to Hiroshi Watanabe, PhD, Computational Biomechanics Div., Institute of Environmental Studies, Graduate School of Frontier Sciences, The University of Tokyo Hongo 7-3-1, Bunkyo-ku Tokyo 113-0033, Japan. Tel.: +81-3-5841-8589; Fax: +81-3-3818-0835; E-mail: nabe{at}sml.k.u-tokyo.ac.jp.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
To relate the subcellular molecular events to organ level physiology in heart, we have developed a three-dimensional finite-element-based simulation program incorporating the cellular mechanisms of excitation-contraction coupling and its propagation, and simulated the fluid-structure interaction involved in the contraction and relaxation of the human left ventricle. The FitzHugh-Nagumo model and four-state model representing the cross-bridge kinetics were adopted for cellular model. Both ventricular wall and blood in the cavity were modeled by finite element mesh. An arbitrary Lagrangian Eulerian finite element method with automatic mesh updating has been formulated for large domain changes, and a strong coupling strategy has been taken. Using electrical analog of pulmonary circulation and left atrium as a preload and the windkessel model as an afterload, dynamics of ventricular filling as well as ejection was simulated. We successfully reproduced the biphasic filling flow consisting of early rapid filling and atrial contraction similar to that reported in clinical observation. Furthermore, fluid-structure analysis enabled us to analyze the wave propagation velocity of filling flow. This simulator can be a powerful tool for establishing a link between molecular abnormality and the clinical disorder at the macroscopic level.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
The heart is an essential organ for maintaining human life by propelling blood throughout the circulatory system. Although recent studies in cellular and molecular biology have greatly promoted our understanding of the structure and function of the heart at the microscopic level, these approaches have contributed little to the understanding of the integrated function of this organ as a pump. In this regard, integration based on computational sciences can be a powerful tool for moving heart research into a new era (Hunter and Borg, 2003Go; Noble, 2002Go).

To construct a comprehensive model of the heart that can accurately simulate the series of events during cardiac cycle, microscopic as well as macroscopic mechanisms should be taken into consideration. Furthermore, description of each event involves the coupling of various disciplines such as electricity, physical chemistry, solid mechanics, and fluid dynamics (multiphysics simulation). Among these, many simulation studies have already been reported in the field of electrophysiology both at the cellular (Luo and Rudy, 1991Go; Noble et al., 1998Go) and organ levels (Kohl et al., 2000Go; Nakazawa et al., 1999Go; Winslow et al., 2000Go) and some of them could successfully integrate the cellular models to simulate the genesis and evolution of arrhythmia in the whole heart (Nakazawa et al., 1999Go; Winslow et al., 2000Go). Similarly in the field of mechanics, attempts have been made to describe the cross-bridge kinetics (Negroni and Lascano, 1996Go; Peterson et al., 1991Go; Rice et al., 1999Go), motion of the left ventricle (LV) muscle (Hunter et al., 1988Go; Huyghe et al., 1992Go), and fluid-structure interaction in the ventricle (Kovacs et al., 2001Go; McQueen and Peskin, 2000Go). However, so far, none of these models have simulated blood flow and the wall motion of the heart powered by the subcellular molecular mechanisms including electrical activity.

To simulate the pumping action of the human LV, we have already developed a fluid-structure interaction (FSI) finite element (FE) model incorporating the propagation of excitation and subcellular excitation-contraction (E-C) coupling mechanisms of individual cardiac myocyte (Watanabe et al., 2002Go). In addition, in this model, an arbitrary Lagragian Eulerian (ALE) finite element method (FEM) with automatic mesh updating has been formulated for large domain changes, and a strong coupling strategy has been taken. As a result, intraventricular blood flow distribution (macroscopic finding) could successfully be related to the molecular mechanism of cardiac contraction for the first time. However, because this was the simple prototype consisting of only LV and systemic arterial system as an afterload, further improvements to the model were necessary to simulate the various aspect of the heart in both normal and diseased states.

Accordingly, in this study, we have incorporated the dynamics of the left atrium (LA) and pulmonary circulation into the model to simulate the ventricular filling dynamics, one of the important issues in cardiology. The results were compared with both experimental and clinical findings to show the usefulness of this method. Especially, information on intraventricular flow provided by the model facilitated the detailed comparison with the observation by Doppler echocardiography including color M-mode Doppler measurement (Brun et al., 1992Go).


   METHOD
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
Details of the formulations adopted for the model have been described in our previous report (Watanabe et al., 2002Go). Here, we briefly present the framework of the model.

Excitation and propagation
Upon stimulations, myocardial cells are excited and contract to develop force. In this model, all of the myocardial cells are electrically connected so that the excitation initiated at a single site can propagate to all of the myocardial cells. By combining the FitzHugh-Nagumo (FHN) model (FitzHugh, 1961Go; Nagumo et al., 1962Go), with the monodomain propagation model (Hodgkin and Huxley, 1952Go), the excitation and its propagation in myocardial tissue can be given by

(1)

(2)
where x and y are variables related to the membrane potential and the adiaphoria, respectively, z is the external stimulation, and g, a, b, and c are constants. The suffix i denotes the spatial position of a unit (a finite element in this article), and j indicates the units surrounding it.

E-C coupling
After excitation (membrane depolarization), a series of subcellular events leads to a transient increase in intracellular calcium concentration [Ca2+], which in turn controls the interaction of contractile proteins namely actin and myosin (cross-bridge kinetics) to develop a force. To describe the dynamic relationship between [Ca2+] and cross-bridge kinetics, a four-state model proposed by Peterson et al. (1991)Go was employed. To connect the membrane depolarization (FHN model) and the four-state model for E-C coupling, a FHN model was used to give a trigger (timing) for the phasic change in Ca2+ ion concentration (Ca2+ transient) described by the following function of time (Peterson et al., 1991Go):

(3)
where [Ca2+]max = 1.2 µM and {tau}Ca = 80 ms (Peterson et al., 1991Go). This equation is analogous to the expression for displacement of a critically damped harmonic oscillator, and thus, this allows for a rapid increase, then a slower decline (Allen and Kurihara, 1982Go; Yue, 1987Go).

The binding of Ca2+ with regulatory protein troponin C (TnC) releases the inhibition to allow the contractile proteins to form "attached" cross-bridges responsible for force generation. However, because it has been shown that the attached cross-bridges remain attached, even after the unbinding of Ca2+ from TnC, the following two states belong to the force-generating populations: 1), TnC binds with Ca2+ and a cross-bridge is attached (O), and 2), TnC does not bind with Ca2+ but the cross-bridge remains attached (O'). In addition, the following two states are defined: 3), TnC binds with Ca2+ but the cross-bridge is detached (U), and 4), TnC does not bind with Ca2+ and the cross-bridge is detached (U'). The complete system can be represented by evolutional equations in matrix form as follows:

(4)
where [Ca2+]f is the concentration of free calcium (not bound to TnC) and k1, k2, k3, k4, ka, kd, and are constants that determine the transition rate between the states. The force is given by the sum of the populations in two "attached" states:

(5)

The calculated time courses of normalized membrane potential xN = x/xmax in the FHN model, [Ca2+]f, and F are shown in Fig. 1 A. In this simulation, F was modified so as to incorporate shortening deactivation (Leach et al., 1980Go) due to the ejection of blood starting at 0.1 s.



View larger version (28K):
[in this window]
[in a new window]
  		FIGURE 1  Characteristics of cell model (A) dynamic relations among normalized membrane potential (dashed line), intracellular calcium concentration [Ca2+] (thin solid line), and normalized developed force (thick solid line) as functions of time. (B) Time-varying stress-strain relationship in fiber direction. As the relative force in A is increased, the stress-strain relation shifts upward left (arrowhead).

 
Constitutive equation of cardiac muscles
To characterize the properties of cardiac muscles, we adopted the Lin-Yin model (Lin and Yin, 1998Go) for the constitutive equation based on hyperelastic material theory. In this model, the strain energy potential W is divided into two components: one is passive (Wpass), and the other is active (Wact):

(6)

(7)

(8)

In the above equations, I1 is the first strain invariant defined as I1 = tr C, where C is the right Cauchy-Green deformation tensor. I4 is a parameter that represents the deformation in muscle fiber direction defined as I4 = nCn, where n is a unit vector defining the preferred direction of muscle fibers in the undeformed state. Under the assumption of incompressibility, the second Piola-Kirchhoff stress is derived as follows:

(9)

(10)

Because the mechanical properties of cardiac muscle vary depending on the population of the attached cross-bridges, F, we made the coefficients for the active components, c5 ~ c10, to be functions of F, whereas those for the passive components, c1 ~ c4, to be constants. To model the phasic change in the stress-strain relationship during the cardiac cycle, the coefficients were expressed as

(11)
where FR is the relative force as defined by FR = F/Fmax. In Fig. 1 B, the calculated stress-strain relationships in uniaxial state are shown for various FR.

FSI analysis
Here, we briefly describe the outline of the fluid-structure interaction FE analysis (for precise formulation; see Zhang and Hisada, 2001Go). Also a short review of ALE method and strong coupling method are shown in Appendix). The fluid is assumed to be Newtonian. The ALE form of the Navier-Stokes equations are discretized using the Galerkin method and FEM. Upon total Lagrangian formulation, the equilibrium equations for structure are discretized using FEM. The full interaction conditions between the fluid and structure are introduced, i.e., the geometrical compatibility conditions and the equilibrium conditions are given. These conditions are automatically satisfied through the element assemblage process. Thus, the system equations in the strong coupling method are constructed. For the time integration method, the extended predictor-multicorrector algorithm is employed based on the Newmark-ß time integration scheme.

In this simulation, the boundaries of the fluid can be assumed to consist of open boundaries and the interface with a deformable structure. Under these boundary conditions, the following methods are employed for mesh control in the fluid domain. The mesh is fixed in space for open boundaries. On the interface with a deformable structure, due to the nature of the strong coupling method, the Lagrangian motion is enforced. In other words, the ALE coordinates are set to be identical to the Lagrangian coordinates on the deformable interface. With these boundary conditions, the motion of the interior nodes is determined by solving the Laplace equation.


   FE MODEL OF THE LV
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
Initial shape of the LV and FE mesh
Although it is an essential requirement for FE analysis to determine the mesh of an object based on the initial shape in an undeformed state under stress-free conditions, it is difficult to know this for a living human LV. Accordingly, we estimated the initial shape based on the following assumptions:

  1. The cavity volume of an adult human LV is 120 ml when the pressure of LV (LVP) is 1500 Pa under fully relaxed conditions (at end diastole).
  2. Because it has been shown that if isovolumic relaxation phase is forced to proceed at the end-systolic volume (ESV) (normally 40–45% of the end-diastolic volume (EDV), thus being ~50–55 ml in this model), the LVP becomes negative to reach the negative part of the end-diastolic pressure volume relationship. We, therefore, assumed the initial volume of LV (LVV) (LVV at which LVP is 0 Pa) is slightly larger than these ESV values, i.e., 60 ml.

To create an FE mesh of the initial shape of the LV (volume 60 ml), we first generated a mesh from the magnetic resonance (MR) images of an adult human LV at end diastole, then modified it by applying an appropriate negative pressure. The resultant FE mesh is shown in Fig. 2.



View larger version (64K):
[in this window]
[in a new window]
  		FIGURE 2  Diagram of macroscopic model. FEM meshed for solid (ventricular wall) and fluid (blood) elements are shown with the electrical analog of afterload (systemic arterial tree) and preload (pulmonary circulation with active LA). PV, pulmonary source pressure; RV, source resistance; CP, pulmonary venous capacitance; RP, pulmonary resistance; LP, pulmonary inertance; RAV, atrioventricular resistance; LAV, atrioventricular inertance; FRL, pulmonary venous flow; FMi, blood flow through the mitral valve; qP and qLA, volumes of the pulmonary venous and LA.

 
It is known that the ventricular wall consists of multiple layers of muscle bundles (Yoran et al., 1973Go). In this analysis, the ventricular wall was divided into six layers from endocardium to epicardium, the fiber directions of which were –60, –30, 0, 30, and 60°, respectively, relative to the plane perpendicular to the long axis of the LV. As for the inflow and outflow of the LV, the following points were taken into consideration. According to the MR imaging data of Kilner et al. (2000)Go, the aorta (Ao) is connected to the LV at an angle of 40° with respect to the long axis of the LV, and the direction of blood flow though the mitral valve is tilted toward the lateral wall by 10° relative to the long axis of the LV. We modeled the inflow and outflow based on these data with a short segment of outflow tract to stabilize the blood flow ejected from the LV (Fig. 2). Also the inflow and outflow velocity profile at the open boundary is assumed constant, i.e., all nodal velocity forced to take the same value.

We used Q1-P0 mixed hexahedral elements (eight nodes for the bilinear velocity or displacement interpolation/constant pressure field) for both the fluid and the structure. To evaluate the spatial convergence of the model, we repeated calculation using two kinds of meshes with different spatial resolution and compared the results. The coarse mesh consisted of 2636 solid and 2448 fluid elements having ~20,000 no. of degree of freedom (NDOF). The fine mesh consisted of 18,976 solid and 9792 fluid elements having ~120,000 NDOF.

In the electrophysiological analysis, the same structural mesh was used as an electrophysiological unit. In other words, each FHN cell unit is embedded in each finite element, and it is electrically connected with the cells in the surrounding elements (see Fig. 2). The excitation of one element, thus, conducts to surrounding elements under the assumption of isotropic conduction.

Models of systemic arterial tree and pulmonary circulation including LA
To simulate the systemic arterial tree, the windkessel model (Westerhof et al., 1971Go) characterized by the following equations was used (see Fig. 2).

(12)

(13)

(14)

P is the aortic pressure, FAo is the blood flow ejected from the LV, R1 is the characteristic impedance, R2 is the peripheral resistance, and C is the capacitance. In this analysis, the following values are employed for these parameters: R1 = 1.00 MPa s/m3, R2 = 150 MPa s/m3, and C = 0.016 µm3/Pa (Murgo et al., 1980Go; Stergiopulos et al., 1999Go).

To model the dynamics of the LA and pulmonary circulation, an electrical analog model proposed by Alexander et al. (1987)Go was used with some modifications (Eqs. 15–18).

(15)

(16)

(17)

(18)

Atrial function is usually characterized in terms of reservoir, conduit, and pump function. In this model, the LA was characterized by the time-varying elastance defined as where PLA is pressure of LA (LAP), qLA is volumes of LA (LAV), and q0 is the volume axis intercept. To model the reservoir and pump function of LA, the time course of ELA (time varying elastance of LA) was defined as follows: 1), in the resting state in which LA acts as a reservoir, LAP was a linear function of LAV, i.e., ELA(t) was constant, ELA(t) = 13.3 MPa/m3; and 2), the contraction of LA started 0.15 s earlier than that of the LV (atrio-ventricular conduction delay) and lasted for 0.2 s with a time course approximated by a sinusoid, MPa/m3, where t0 is cardiac cycle. PV is the pulmonary source pressure, RV is the source resistance, CP is the pulmonary venous capacitance, RP is the pulmonary resistance, LP is the pulmonary inertance, RAV is the atrioventricular resistance, LAV is the atrioventricular inertance, FRL is the pulmonary venous flow, FMi is the blood flow through the mitral valve, qP is the volume of the pulmonary venous, and PLV is the pressure boundary condition for the FSI model (see Fig. 2). The following values are assumed for these parameters: PV = 1600 Pa, RV = 0.333 MPa s/m3, CP = 0.0525 µm3/Pa, RP = 12.0 MPa s/m3, LP = 0.133 MPa s2/m3, RAV = 3.47 MPa s/m3, LAV = 0.267 MPa s2/m3.

Material properties
The coefficients in the passive part of the Lin-Yin model (c1–c4) were determined so that the LVV becomes 120 ml when the LVP is 1500 Pa under fully relaxed conditions. This was achieved by using 80% of the average value reported by Lin and Yin (Lin and Yin, 1998Go). The coefficients in the active part (c5–c10) were also determined to meet the following conditions:

  1. Isolated cardiac myocytes contract up to 20% under an unloaded condition.
  2. The maximum value of the time derivative of the LVP (max dP/dt) is between 200 kPa/s (1500 mmHg/s) and 240 kPa/s (1800 mmHg/s) in isovolumic contraction phase.

First, we determined the c10 required to fulfill condition 1. Then, c5 ~ c9 were determined to achieve condition 2. The coefficients thus determined were as follows: c1 = 82.4 Pa, c2 = 7.30, c3 = 1.86, c4 = 0.0640, In addition, the density of cardiac muscle was assumed to be {rho} = 1.37 x 103 kg/m3, the density of blood {rho} = 1.06 x 103 kg/m3, and the viscosity of blood µ = 4.71 x 10–3 Pa s.

Simulation
In the case of normal contraction, the excitation conducting from the atria spreads rapidly over the entire endocardium via the conduction system (Purkinje fibers) and then propagates toward the epicardium within a short period of time (<40 ms). To meet this condition, Eqs. 1 and 2 are solved with the following nondimensional initial values and coefficients: x = 1.20, y = 0.624, a = 0.70, b = 0.80, c = 0.70, g = 1.00, and z = 20.0, and the stimulation time of 50 ms was imposed on all of the elements in the endocardium.

The following boundary and initial conditions were set for each period during the cardiac cycle. During the isovolumic contraction period, both the mitral and aortic valves were closed so that the flow velocity was zero at those two boundaries. The initial value of LVV was 120 ml at an LVP of 1500 Pa under fully relaxed conditions. Then, the muscles were activated by the excitation-contraction mechanism (represented by the FHN model and the four-state model) to increase the LVP. When the LVP exceeded the aortic pressure (AP) calculated by the windkessel model, the aortic valve opened. During the ejection period, with the mitral valve closed, the pressure boundary condition given by the windkessel model was applied to the open aortic valve area. The FSI FE model and the windkessel model were iteratively calculated to determine the equilibrium state, and when backward flow occurred in the aorta, the aortic valve closed to finish ejection. From this point, the windkessel was disconnected from the LV and its status was calculated under the condition that FAo = 0 until the next cardiac cycle. During isovolumic relaxation period, the relaxation of cardiac muscles proceeded while both the aortic and mitral valves closed. When the LVP fell below the LAP, the mitral valve opened, and the filling period began. Based on these assumptions, we simulated three cycles of contractions for a heart rate (HR) of 60, 75, 90, and 120 beats per minute (one cardiac cycle: 1.0 s, 0.8 s, 0.665 s, and 0.5 s, respectively). In this model, Ao and LA are not modeled by finite element, so that the displacement and the velocity of the node located in the atrio-ventricular ring are fixed to zero though the cardiac cycle. The time step is set to 1.0 ms during the ejection phase, and 2.5 ms in other phases.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
Convergence of the model
Aortic and mitral flows during the cardiac cycle calculated with coarse and fine meshes are compared in Fig. 3 A. Traces for two kinds of meshes were almost superimposable with the maximum aortic flow being 464 ml/s (coarse) and 484 ml/s (fine), and the maximum mitral flow 192 ml/s (coarse) and 195 ml/s (fine). Similarly, LV, aortic, and LA pressures agreed reasonably well between the two kinds of meshes (Fig. 3 B). From these results, it is concluded that the results with the coarse mesh had enough accuracy but we presented the results obtained from the fine mesh in the following text so as to understand the flow field visually better. However, because, we consider that the current mesh discretization is still insufficient for resolving the boundary layer of fluid, we are trying to further improve the model with the higher performance computers and parallel processing.



View larger version (16K):
[in this window]
[in a new window]
  		FIGURE 3  Comparison of two kinds of meshes with different spatial resolution (A) LV inflow and outflow are shown for the fine (solid line) and the coarse (dashed line) meshes. (B) LVP (thick lines), AP (hair lines), and LAP (thin lines) are shown for the fine (solid lines) and the coarse (dashed lines) meshes.

 
Activation and intraventricular flow pattern
Time-lapse images of LV contraction and resultant blood flow are shown in Fig. 4 together with the activation of myocardial tissue (strain in the fiber direction). (Also see Supplementary Material.) Myocardial strain develops initially at the endocardium then spreads toward the epicardium and relaxes from the endocardium (0.4 s). Compared to the slow propagation in strain development, relaxation was relatively fast, because, in this simulation, relaxation was accelerated by the shortening deactivation. It is of note that blood flow continues even in isovolumic period (both valves closed) and forms a large vortex (0.0–0.1 s). This phenomenon has been recently reported by MR imaging technique and considered to facilitate blood ejection during the following period (Kilner et al., 2000Go). To facilitate the understanding on flow dynamics in the ventricle, LVP, and flow distributions are shown at the beginning of contraction (Fig. 5, left), during ejection (Fig. 5, middle), and during filling period (Fig. 5, right). The flow continues even where the pressure gradient is minimal indicating that the intraventricular flow pattern is governed also by the momentum of the blood. Those figures also confirm the fact that the pressure is low where fast flow is observed (Bernoulli's law).



View larger version (79K):
[in this window]
[in a new window]
  		FIGURE 4  Time-lapse images of LV contraction and intraventricular blood flow. Numbers in each plate indicate the time after the onset of activation in seconds. Activation level (relative force) is shown in color coding. Velocity of blood is shown as the length of the arrows.

 


View larger version (28K):
[in this window]
[in a new window]
  		FIGURE 5  Intraventricular flow and pressure distribution. Numbers in each plate indicate the time after the onset of activation in seconds. Pressure level is shown in color coding. Velocity of blood is shown as the length of the arrows.

 
Wall motion
Changes in geometry including the wall motion during the LV contraction and relaxation in Fig. 4 are summarized in Fig. 6. In clinical cardiology, atrio-ventricular valve plane displacement and its velocity are used as important indices of pump function of LV (Alam and Hoglund, 1992Go; Hoglund et al., 1988Go; Lundback, 1986Go). However, because this model does not include LA and/or Ao as structural entities, the LV was fixed in space at the atrio-ventricular valve plane. Therefore, we analyzed the long axis displacement by measuring the vertical displacement of the apex as an alternative (Fig. 6 A). As reported using echocardiography, the model ventricle shortened during ejection and lengthened during filling, but the maximum long axis displacement of LV (~10 mm) was definitely smaller compared to the reported value for normal subjects (~16 mm) (Alam and Hoglund, 1992Go; Hoglund et al., 1988Go; Lundback, 1986Go). If we look closer at the time course (Fig. 6 A), we can identify a hump immediately after the lengthening caused by the rapid filling (0.8–1.0 s). This was caused by the rebound shortening of the model ventricle without damping property and could be the cause of smaller vertical displacement. The lack of damping property also caused oscillations of relatively high frequency (~20 Hz). during ejection (Fig. 6, A and B). Ventricular torsion is shown in Fig. 6 B. During contraction, LV initially twisted clockwise by a small amount, then a counter-clockwise direction as it ejected the blood. Moon et al. studied the torsion of canine LV using radio-opaque marker technique and reported the similar time course of torsional deformation with the maximum amplitude of ~6° at baseline and ~17° under inotropic intervention (Moon et al., 1996Go). This result also agreed with this observation. Wall thickness and its time derivative (rate of thickness change) are shown in Fig. 6 C. During systole, wall thickness increased by ~5mm. It is well known that, at end systole, the wall thickness of septum is twice as large as that of posterior wall (Lundback, 1986Go). The value obtained from the current model with axisymmetric structure was between those reported for septum and posterior wall. During the filling period, wall thickness rapidly decreased with the peak rate of ~3.2 cm/s. This value is smaller than that of diastolic posterior wall velocity for normal subject by Fujii et al. (7.3 cm/s) (Fujii et al., 1979Go) but definitely larger than that for interventricular septum estimated from the Fig. 1 of the same report. We consider that this was also due to the axisymmetric shape of the model ventricle.



View larger version (16K):
[in this window]
[in a new window]
  		FIGURE 6  Wall motion of the simulated LV. (A) Vertical displacement, (B) rotational movement, (C) wall thickness (solid line), and changing rate of wall thickness rate (dashed line).

 
Hemodynamic parameters and filling dynamics
Hemodynamic parameters during the LV contraction and relaxation in Fig. 3 are summarized in Fig. 7. As is shown in the pressure-volume diagram (Fig. 7 A), the simulated LV achieved normal ejection fraction (55%) with 66 ml of stroke volume (corresponding cardiac output = 0.066 x 60 = 4.0 L/min). Nonlinear pressure-volume relation of LV during diastole was also reproduced. Fig. 4 B shows the changes in LVP, LAP, and AP during cardiac cycle. It is of note that, during ejection period, the LVP was higher than the AP at the beginning but the AP exceeded the LVP in the latter part of the phase, although ejection still continued (Fig. 7 B, inset). This phenomenon has been reported in animal experiments and considered to be caused by the inertia of the blood (Noble, 1968Go). At the time of mitral valve opening, LA pressure is higher than LV pressure but they gradually equilibrate as the LV is filled with the blood from LA and LV pressure exceeds in mid-diastole. At end-diastole, LA pressure again becomes higher because of the active contraction of LA. These changes in trans-mitral pressure differences accelerate and decelerate the blood flow across mitral valve to result in early rapid filling (E) wave and atrial contraction (A) wave shown in Fig. 7 C. Obviously, the amplitude of E-wave is greater than A-wave as is the case with normal heart. All these results are in good agreement with experimental studies (Ishida et al., 1986Go; Steine et al., 2002Go). The change in LVV induced by these outflow and inflow is shown with the concomitant change in LAV (Fig. 7 D). The volumes of these two chambers changed in a symmetrical manner but the sum of them did not remain constant because of the direct flow from pulmonary circulation to LV (conduit function of LA). When normalized to the end diastolic value, the sum (LVV + LAV) varied by ~25% with the maximum deviation occurring at end systole, which agreed well with reported data for normal subjects (Bowman and Kovacs, 2004Go).



View larger version (22K):
[in this window]
[in a new window]
  		FIGURE 7  Hemodynamics of the simulated LV. (A) Pressure-volume diagram, (B) dynamic changes in LVP (thin solid line), AP (dashed line), and LAP (thick solid line). (Insets) Magnified views of the relation between (1) LVP and AP during ejection (2) LVP and LAP during filling. (C) Aortic flow (thick solid line) and mitral flow (thin solid line). (D) Time courses of LVV (thin solid line), LAV (thick solid line), and total volume of LV and LA (dashed line).

 
The FSI analysis could simulate another important index of LV diastolic filling. Recently, it is proposed that the diastolic flow wave front velocity measured by color M-mode echocardiography along the line connecting the mitral valve and apex (wave propagation velocity) can be a useful index of ventricular relaxation property. According to the practice, we plotted the color-coded flow velocity data at each sampling point (measured as the distance from the mitral valve) on the y axis as a function of time in Fig. 8. The flow wave front during early filling period is identified as a linear margin and the flow propagation velocity is defined as the slope of this linear segment. The wave propagation velocity thus evaluated is ~0.8m/s, which again agreed well with reported data for normal subjects (Brun et al., 1992Go).



View larger version (60K):
[in this window]
[in a new window]
  		FIGURE 8  Color-coded map of the flow velocity along the long axis of the LV. At each time point (ordinate), flow velocity is shown in color at each distance from the mitral valve (abscissa). The slope of the thick line connecting wave front is the wave propagation velocity. Velocity level is shown in color coding.

 
Effect of heart rate
To study the effect of heart rate, we performed the simulation at 60, 75, 90, and 120 beats/min and compared the results in Fig. 9. As the heart rate increased, E- and A-waves of diastolic filling merged and became indistinguishable whereas the aortic flow remained relatively constant until the heart rate reached 120. These results are qualitatively similar to clinical observations but the transition occurred at a much lower heart rate. Possible reasons for this discrepancy are: 1), durations of contraction of LV and LA was fixed, and 2), atrio-ventricular conduction delay was fixed at 0.15 s. In physiological settings, both of these parameters were controlled by neuronal and hormonal regulations to compensate for the curtailed time for diastolic filling.



View larger version (16K):
[in this window]
[in a new window]
  		FIGURE 9  Aortic and mitral flows at different HR (thin solid line, HR = 60 bpm; thin dashed line, HR = 75 bpm; thick solid line, HR = 90 bpm, thick dashed line, HR = 120 bpm).

 
Comparison with other methods
Detailed models of the ventricles having fiber orientations and sheet structure have already been constructed (Costa et al., 2001Go; LeGrice et al., 2001Go; Nash and Hunter, 2000Go). By coupling these models with cellular models of electrophysiology, propagation of activation was simulated for normal and ectopic excitation (Kohl et al., 2000Go; Winslow et al., 2000Go). These anatomical models were also used to establish relations between global hemodynamics to regional mechanics of the myocardial tissue (Nash and Hunter, 2000Go; Stevens and Hunter, 2003Go; Usyk et al., 2000Go). Although these researchers successfully reproduced the three-dimensional strains (Nash and Hunter, 2000Go; Usyk et al., 2000Go) and sarcomere length (Stevens and Hunter, 2003Go) at various sites in the ventricular wall during both diastole and systole, these simulations differ from this model in the following respects: 1), tissue properties are characterized by constitutive equation that bears no relation to the subcellular mechanisms (Costa et al., 2001Go; Huyghe et al., 1992Go) or, even if related, which are not dynamic functions of time (only at end-systole) (Nash and Hunter, 2000Go); and 2), intraventricular flow and its interaction with the ventricular wall were not take into consideration. As a result, these models do not provide us with any information on the flow dynamics during ventricular ejection and filling. Although Peskin and colleagues (Kovacs et al., 2001Go; McQueen and Peskin, 2000Go), using immersed boundary method, reported the detailed flow mapping in the four chambers of the heart and great vessels, their methodology, in principle, does not allow one to incorporate the subcellular mechanisms to the material property.

Simulation of diastolic filling
With the advancement of Doppler echocardiography, a pattern of LV filling has been analyzed to detect the diastolic dysfunction. Although, so far, several indices have been proposed such as the peak velocity of early rapid filling, atrial contraction, their ratio, and the acceleration and deceleration times of the early filling wave, they are mostly empirical in nature (Thomas and Weyman, 1991Go). Accordingly, simulation (Kovacs et al., 1987Go; Meisner et al., 1991Go; Thomas and Weyman, 1991Go, 1992Go) as well as experimental studies (Ishida et al., 1986Go) has been performed to establish a link between the tissue property and filling pattern. Thomas and Weyman (1991Go, 1992Go), using a lumped parameter model, successfully simulated the influence of the left ventricular relaxation rate on the peak velocity of early rapid filling, atrial contraction, and their ratio, but it is beyond the scope of the lumped model to examine the effect of the propagation of excitation and/or detailed cellular function. Furthermore, flow propagation velocity, which is now emerging as a useful index of ventricular relaxation rate (Brun et al., 1992Go), can only be simulated by the three-dimensional fluid-structure-interaction model in this study. We are now attempting to simulate the filling patterns in various diseased states including cardiomyopathy and ischemia.

Study limitation
Toward the complete modeling of the left ventricular filling, there are still many aspects to be added to the current model. Because of the limitation in computational power, LA was modeled as an electrical analog circuit. This treatment might have obscured the physiology of LA, which is usually characterized by reservoir, conduit, and pump function. However, from the analogy to LV time-varying elastance model (Suga and Sagawa, 1974Go), pump function can be characterized by the maximum value of elastance of LA (ELA(t) in Fig. 2). Reservoir function is evaluated by the difference between the maximum and minimum LA volumes (LA reservoir volume), which is influenced by the minimum value of ELA(t). Conduit function can also be evaluated by the difference between LV stroke volume and LA reservoir volume in Fig. 7 D. However, to study the ventricular stroke under a variety of conditions ranging from normal to abnormal, LA with realistic shape and tissue property is necessary.

The constitutive equation proposed by Lin and Yin (1998)Go was employed to model the myocardial tissue in which material property was divided into passive and active parts. Active part represents the change in myocyte property induced by the cross-bridge formation. Passive property of the myocardium is determined by the extracellular matrix (collagen) and passive components of myocyte (cytoskeleton including titin). Considering the complexity and importance of passive property of these components (Hein et al., 2002Go), a detailed model is required to simulate various diseased state such as cardiac hypertrophy; however, in this model, all these parameters are lumped in the anisotropic stress-strain relations.

It is recognized that diastolic function is determined by dynamic and passive chamber stiffness and relaxation of the ventricle. In this model, relaxation is characterized by the time constant for the decay of Ca2+ transient ( in Eq. 3). Passive chamber stiffness is formulated by the constitutive equation representing the passive component (Eq. 6). However, we did not incorporate the dynamic property of chamber stiffness in the model. Although it has been shown that the difference between dynamic and passive diastolic property is slight in a normal subject, there is a marked discordance in diseased states such as hypertrophic cardiomyopathy (Pak et al., 1996Go). Dynamic property should also be incorporated in the future modeling.

The FitzHugh-Nagumo model was adopted for simulating electrical activity in this model. As this is obviously an oversimplification of reality in which many transmembrane ionic currents (Na, K, Ca, Cl, etc.) and intracellular [Ca2+] dynamics determine depolarization, we are now improving the model by incorporating ion channel models proposed by Faber and Rudy (2000)Go with 16 membrane currents.


   SUMMARY
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
To simulate the contraction and the filling dynamics of LV, an ALE fluid-structure interaction FE code was developed and coupled with the electrical analog of the LA and pulmonary circulation. Although the descriptions of each component were not novel, to our knowledge, this is the first report in which subcellular molecular events were integrated to produce macroscopic blood flow that routine clinical tests can reveal for each individual. The results of the simulation reproduced well the clinical observation concerning the indices of diastolic function. With further improvements, the simulator can be truly a practical tool for the study of heart disease and clinical diagnosis.


   APPENDIX
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
This section briefly describes the ALE method and strong coupling method. (for precise formulation, see Zhang and Hisada, 2001Go). Based on the work by Huerta and Liu (1988)Go, the ALE method is formulated for a deformable fluid domain as follows. We assume the viscous fluid to be isothermal and that, with B, P, and {rho} being the fluid bulk modulus, the pressure, and the fluid density, respectively. The ALE form of the Navier-Stokes equations can be expressed as

(A1)

(A2)
where is the time derivative in the ALE coordinates, the velocity of which is vm, and the convective velocity c = v – vm. {sigma} is the Cauchy stress and g the acceleration of gravity. denotes the spatial domain with the boundary at time t. The superscript f stands for the fluid component. The fluid is assumed to be Newtonian. Using the Galerkin method and finite element discretization for Eqs. A1 and A2, we obtain the following equations in the matrix form:

(A3)

(A4)
where and M are the generalized mass matrices for pressure and velocity, respectively, and and {Lambda} are the generalized matrices of convective terms for pressure and velocity, respectively. is the fluid viscosity matrix, G is the divergence operator matrix, and are the time derivatives of the pressure and velocity vectors, respectively, in the ALE coordinates, and F is the external force vector. As for the stabilization, standard streamline upwind/Petrov Galerkin formulation (Tezduyar et al., 1992Go) is employed.

Let be the spatial domain of a structure. Here the superscript s stands for the structural component. The equilibrium equations for structure are:

(A5)

Upon total Lagrangian formulation and finite element discretization, a nonlinear system of equations is obtained in the matrix form, and the incremental form at each time step is given as follows:

(A6)

Here, M is the mass matrix, K the tangent stiffness matrix, and F the external force vector. is composed of the internal force and the inertial force of the structure, and and are the increments of the acceleration and the displacement vectors, respectively.

The full interaction conditions between the fluid and structure are introduced, i.e., the geometrical compatibility conditions:

(A7)
and the equilibrium conditions:

(A8)
are given, where nf and ns are the outward normals corresponding to the fluid and the structure, respectively, and is the interface between the fluid and the structure.

The above conditions are automatically satisfied through the element assemblage process. Thus, the system equations in the strong coupling method are constructed, and the incremental form can be expressed as

(A9)
where and are the variable vectors of the coupled system, as defined by

(A10)
Here, P is the pressure of the fluid, the velocity vector of the fluid independent of the structure, the coupled velocity vector, the velocity vector of the structure independent of the fluid, the coupled displacement vector, and the displacement vector of the structure independent of the fluid. In Eq. A9, F denotes the external force vector of the coupled system and the equivalent internal force vector including all effects. denotes the mass matrix composed of those for the fluid and the structure. consists of the divergence, viscous, and convective terms of the fluid. is the tangent stiffness matrix of the structure.


   SUPPLEMENTARY MATERIAL
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
ACKNOWLEDGEMENTS
 REFERENCES
 
An online supplement to this article can be found by visiting BJ Online at http://www.biophysj.org.


   ACKNOWLEDGEMENTS
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
REFERENCES
 
The authors thank Dr. K. Sunagawa, Dr. M. Sugimachi, Dr. K. Nakazawa, and Dr. M. Inagaki at the National Cardiovascular Research Center for their valuable comments.

This work was supported in part by the grants from Core Research for Evolutional Science and Technology of the Japan Science and Technology Agency and National Cardiovascular Research Center, The Program for Promotion of Fundamental Studies in Health Sciences of the Organization for Pharmaceutical Safety and Research, Suzuken Memorial Foundation, and The Research Grant for Cardiovascular Disease from the Ministry of Health, Labor and Welfare.

Submitted on October 16, 2003; accepted for publication May 28, 2004.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 METHOD
 FE MODEL OF THE...
 RESULTS AND DISCUSSION
 SUMMARY
 APPENDIX
 SUPPLEMENTARY MATERIAL
 ACKNOWLEDGEMENTS
 REFERENCES
 
Alam, M., and C. Hoglund. 1992. Assessment by echocardiogram of left ventricular diastolic function in healthy subjects using the atrioventricular plane displacement. Am. J. Cardiol. 69:565–568.[Medline]

Alexander, J., K. Sunagawa, N. Chang, and K. Sagawa. 1987. Instantaneous pressure-volume relation of the ejecting canine left atrium. Circ. Res. 61:209–219.[Abstract/Free Full Text]

Allen, D. G., and S. Kurihara. 1982. The effects of muscle length on intracellular calcium transients in mammalian cardiac muscle. J. Physiol. 327:79–94.[Abstract/Free Full Text]

Bowman, A. W., and S. J. Kovacs. 2004. Left atrial conduit volume is generated by deviation from the constant-volume state of the left heart: a combined MRI-echocardiographic study. Am. J. Physiol. 286:H2416–H2424.

Brun, P., C. Tribouilloy, A. M. Duval, L. Iserin, A. Meguira, G. Pelle, and J. L. Dubois-Rande. 1992. Left ventricular flow propagation during early filling is related to wall relaxation: a color M-mode Doppler analysis. J. Am. Coll. Cardiol. 20:420–432.[Abstract]

Costa, K. D., J. W. Holmes, and A. D. McCulloch. 2001. Modelling cardiac mechanical properties in three dimensions. Philos. Trans. R. Soc. Lond. A. 359:1233–1250.[CrossRef]

Faber, G. M., and Y. Rudy. 2000. Action potential and contractility changes in [Na+]i overloaded cardiac myocytes: a simulation study. Biophys. J. 78:2392–2404.[Abstract/Free Full Text]

FitzHugh, R. 1961. Impulses and physiological states in theoretical models of nerve membrane. Biophys. J. 1:445–466.[Abstract/Free Full Text]

Fujii, J., H. Watanabe, S. Koyama, and K. Kato. 1979. Echocardiographic study on diastolic posterior wall movement and left ventricular filling by disease catogory. Am. Heart J. 98:144–152.[CrossRef][Medline]

Hein, S., W. H. Gaasch, and J. Schaper. 2002. Giant molecule titin and myocardial stiffness. Circulation. 106:1302–1304.[Free Full Text]

Hodgkin, A. L., and A. F. Huxley. 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117:500–544.[Free Full Text]

Hoglund, C., M. Alam, and C. Thorstrand. 1988. Atrioventricular valve plane displacement in healthy persons. Acta Med. Scand. 224:557–562.[Medline]

Huerta, A., and W. K. Liu. 1988. Viscous flow with large surface motion. Comput. Methods Appl. Mech. Engrg. 69:277–324.[CrossRef]

Hunter, P. J., and T. K. Borg. 2003. Integration from proteins to organs: the physiome project. Nat. Rev. Mol. Cell Biol. 4:237–243.[CrossRef][Medline]

Hunter, P. J., A. D. McCulloch, P. M. F. Nielsen, and B. H. Smaill. 1988. A finite element model of passive ventricular mechanics. In Computational Methods in Bioengineering. R. L. Spilker and B. R. Simon, editors. ASME BED, New York, NY. 387–397.

Huyghe, J. M., T. Arts, D. H. Van Champen, and R. S. Reneman. 1992. Porous medium finite element model of the beating left ventricle. Am. J. Physiol. 262:H1256–H1267.[Medline]

Ishida, Y., J. S. Meisner, K. Tsujioka, J. I. Gallo, C. Yoran, R. W. M. Frater, and E. L. Yellin. 1986. Left ventricular filling dynamics: influence of left ventricular relaxation and left atrial pressure. Circulation. 74:187–196.[Abstract/Free Full Text]

Kilner, P. J., G. Z. Yang, A. J. Wilkers, R. H. Mohiaddin, D. N. Firmin, and M. H. Yacoub. 2000. Asymmetric redirection of flow thorough the heart. Nature. 404:759–761.[CrossRef][Medline]

Kohl, P., D. Noble, R. L. Winslow, and P. J. Hunter. 2000. Compuational modelling of biological systems: tools and visions. Philos. Trans. R. Soc. Lond. A. 358:579–610.[CrossRef]

Kovacs, S. J., B. Barzilai, and J. E. Perez. 1987. Evaluation of diastolic function with Doppler echocardiography: the PDF formalism. Am. J. Physiol. 252:H178–H187.[Medline]

Kovacs, S. J., D. M. McQueen, and C. S. Peskin. 2001. Modeling cardiac fluid dynamics and diastolic function. Philos. Trans. R. Soc. Lond. A. 359:1299–1314.[CrossRef]

Leach, J. K., A. J. Brady, B. J. Skipper, and D. L. Millis. 1980. Effects of active shortening on tension development of rabbit papillary muscle. Am. J. Physiol. 238:H8–H13.[Medline]

LeGrice, I., P. J. Hunter, A. Young, and B. Smaill. 2001. The architecture of the heart: a data-based model. Philos. Trans. R. Soc. Lond. A. 359:1217–1232.[CrossRef]

Lin, D. H. S., and F. C. P. Yin. 1998. A multiaxial constitutive law for mammalian left ventricular myocardium in steady-state barium contracture or tetanus. J. Biomech. Eng. 120:504–517.[Medline]

Lundback, S. 1986. Cardiac pumping and function of the ventricular septum. Acta. Physiol. Scand. Suppl. 550:1–101.[Medline]

Luo, C., and Y. Rudy. 1991. A model of the ventricular cardiac action potential depolarization, repolarization, and their interaction. Circ. Res. 68:1501–1526.[Abstract/Free Full Text]

McQueen, D. M., and C. S. Peskin. 2000. A three-dimensional computer model of the human heart for studying cardiac fluid dynamics. Comput. Graph. 34:57–60.[CrossRef]

Meisner, J. S., G. Keren, O. E. Pajaro, A. Mani, J. A. Strom, R. W. Frater, S. Laniado, and E. L. Yellin. 1991. Atrial contribution to ventricular filling in mitral stenosis. Circulation. 84:1469–1480.[Abstract/Free Full Text]

Moon, M. R., A. DeAnda, Jr., G. T. Daughters, N. B. Ingels, and C. Miller. 1996. Effects of chordal disruption on regional left ventricular torsional deformation. Circulation. 94:143–151.[Abstract/Free Full Text]

Murgo, J. P., N. Westerhof, J. P. Giolma, and S. A. Altobelli. 1980. Aortic input impedance in normal man: relationship to pressure wave forms. Circulation. 62:105–116.[Free Full Text]

Nagumo, J., S. Arimoto, and S. Yoshizawa. 1962. An active pulse transmission line simulating nerve axon. Proc. IRE. 50:2061–2070.[CrossRef]

Nakazawa, K., T. Namba, and R. Suzuki. 1999. Dynamics of spiral waves in three dimensional FHN model media: a fundamental study on the mechanism of functional reentrant tachyarrhythmia by computer simulations. Biomed. Eng. 37:63–77.

Nash, M. P., and P. J. Hunter. 2000. Computational mechanics of the heart. J. Elast. 61:113–141.[CrossRef]

Negroni, J. A., and E. C. Lascano. 1996. A cardiac muscle model relating sarcomere dynamics to calcium kinetics. J. Mol. Cell. Cardiol. 28:915–929.[CrossRef][Medline]

Noble, D. 2002. Modeling the heart: from genes to cells to the whole organ. Science. 295:1678–1682.[Abstract/Free Full Text]

Noble, M. I. M. 1968. The contribution of blood momentum to left ventricular ejection in the dog. Circ. Res. 23:663–670.[Abstract/Free Full Text]

Noble, D., A. Varghese, P. Kohl, and P. Noble. 1998. Improved guinea-pig ventricular cell model incorporating a diadic space, IKr and IKs, and length- and tension-dependent processes. Can. J. Cardiol. 14:123–134.[Medline]

Pak, P. H., W. L. Maughan, K. L. Baughman, and D. A. Kass. 1996. Marked discordance between dynamic and passive diastolic pressure-volume relations in idiopathic hypertrophic cardiomyopathy. Circulation. 94:52–60.[Abstract/Free Full Text]

Peterson, J. N., W. C. Hunter, and M. R. Berman. 1991. Estimated time course of Ca2+ bound to troponin C during relaxation in isolated cardiac muscle. Am. J. Physiol. 260:H1013–H1024.[Medline]

Rice, J. J., R. L. Winslow, and W. C. Hunter. 1999. Comparison of putative cooperative mechanisms in cardiac muscle: length dependence and dynamic responses. Am. J. Physiol. 276:H1734–H1754.[Medline]

Steine, K., M. Stugaard, and O. A. Smiseth. 2002. Mechanisms of diastolic intraventricular regional pressure differences and flow in the inflow and outflow tracts. J. Am. Coll. Cardiol. 40:983–990.