Original Research ARTICLE
Toward an integrative computational model of the guinea pig cardiac myocyte
- Department of Biomedical Engineering, Institute for Computational Medicine, The Johns Hopkins University School of Medicine and Whiting School of Engineering, Baltimore, MD, USA
The local control theory of excitation-contraction (EC) coupling asserts that regulation of calcium (Ca2+) release occurs at the nanodomain level, where openings of single L-type Ca2+ channels (LCCs) trigger openings of small clusters of ryanodine receptors (RyRs) co-localized within the dyad. A consequence of local control is that the whole-cell Ca2+ transient is a smooth continuous function of influx of Ca2+ through LCCs. While this so-called graded release property has been known for some time, its functional importance to the integrated behavior of the cardiac ventricular myocyte has not been fully appreciated. We previously formulated a biophysically based model, in which LCCs and RyRs interact via a coarse-grained representation of the dyadic space. The model captures key features of local control using a low-dimensional system of ordinary differential equations. Voltage-dependent gain and graded Ca2+ release are emergent properties of this model by virtue of the fact that model formulation is closely based on the sub-cellular basis of local control. In this current work, we have incorporated this graded release model into a prior model of guinea pig ventricular myocyte electrophysiology, metabolism, and isometric force production. The resulting integrative model predicts the experimentally observed causal relationship between action potential (AP) shape and timing of Ca2+ and force transients, a relationship that is not explained by models lacking the graded release property. Model results suggest that even relatively subtle changes in AP morphology that may result, for example, from remodeling of membrane transporter expression in disease or spatial variation in cell properties, may have major impact on the temporal waveform of Ca2+ transients, thus influencing tissue level electromechanical function.
Since publication of the first computational model of the cardiac myocyte action potential (AP) in 1960 (Noble, 1960), the range of biological processes described in models of the cardiac myocyte has grown continuously. While the integrative nature of today’s most commonly used models differ, the sub-cellular processes for which there are quantitative, experimentally based models include: (a) voltage-gated ion channels and currents; (b) intracellular calcium (Ca2+) dynamics and Ca2+-induced Ca2+-release (CICR); (c) electrogenic and ATP-dependent membrane transporters; (d) regulation of intracellular Ca2+, sodium (Na+), potassium (K+), and hydrogen ion (H+) concentrations; (e) mitochondrial ATP production and its regulation; (f) coupling of ATP production to energy requiring membrane transporters and myofilaments; and (g) ligand gated membrane receptors and intracellular signaling pathways. However few models combine electrophysiology with contraction mechanics, mitochondrial energetics, or intracellular signaling due to the computational difficulty of combining disparate time and/or spatial scales. While the incorporation of more cellular component models increases the descriptive power of the combined model, the complexity of whole-cell model behavior and computational cost increase almost exponentially with the number of constitutive mechanisms represented.
The close interplay between modeling and experiments has enabled a remarkably deep understanding of the function of cardiac myocytes. In some cases, myocyte models not only reconstruct the experimental data on which they are based, they predict new emergent behaviors that have been validated subsequently by experiments. As a result, models now play a central role in understanding the relationships between molecular function and the integrated behavior of the cardiac myocyte in health and disease.
The most fundamental property of cardiac myocytes is that they are electrically excitable cells (they generate APs). The rapid increase in membrane depolarization during the early phase of the AP increases the open probability of sarcolemmal L-Type Ca2+ channels (LCCs). LCCs are preferentially located in the t-tubules (see Figure 1A), which are invaginations of the sarcolemma extending deep into the cell. Further, LCCs within the t-tubules are preferentially localized in structures known as dyads. Dyads are regions where t-tubule membrane is in close opposition to the sarcoplasmic reticulum (SR) membrane. The SR is a luminal organelle located throughout the interior of the cell. It is involved in uptake, sequestration, and release of Ca2+ in a process known as intracellular Ca2+ cycling. The junctional SR (JSR) is the portion of the SR most closely approximating the t-tubules. Ryanodine receptors (RyRs) are channels located in the JSR membrane in close opposition to LCCs in the dyad. During the initial phase of the AP, sarcolemmal membrane depolarization increases LCC open probability. The resulting flux of Ca2+ into the dyadic space (trigger flux) leads to Ca2+ binding to the RyRs. This increases RyR open probability, and when open, the resulting Ca2+ flux through RyRs (release flux) is directed into the dyadic space. There are thousands of dyads within the cardiac myocyte, and the net flux of Ca2+ from dyads into the cytosol triggers muscle contraction in a process known as excitation-contraction coupling (ECC). The ratio of release to trigger flux is typically large, and is referred to as ECC gain.
Figure 1. Overview of guinea pig coupled model. (A) Schematic illustration of the model structure. (B) LCC (top) and RyR (bottom) Markov model structures. (C) The 40-state LCC:RyR model representing each possible pairing of LCC and RyR states.
Graded release refers to the phenomenon, originally observed by Fabiato (1985), whereby Ca2+ release from JSR is a graded, smooth, continuous function of the amount of trigger Ca2+ entering the cell via LCCs. The majority of cardiac myocyte models lump all dyadic spaces together into a common pool known as the subspace, and net trigger flux through LCCs and net release flux through RyRs is directed into this common pool. In a landmark 1992 paper (Stern, 1992), Stern showed that the strong positive feedback effect on RyR open probability due to the fact that release flux is directed into the same pool of Ca2+ that serves as the trigger for RyRs (i.e., the common pool) results in all-or-none rather than graded release. More specifically, he showed that common pool models cannot reproduce both high gain and graded release. Despite this fact, early common pool models were able to reproduce a broad range of cardiac myocyte behaviors. However, Greenstein and Winslow (2002) showed that incorporation of new experimental data, demonstrating that Ca2+-dependent inactivation (CDI) of LCCs is stronger than voltage-dependent inactivation (VDI), into common pool models de-stabilized repolarization of the AP due to abnormal Ca2+ handling (here, de-stabilization means that at normal physiological pacing rates, APs could exhibit an oscillatory plateau phase, and large, irregular variation in AP duration, APD). Thus, common pool models not only fail to capture the graded release property, they cannot capture one of the most fundamental properties of normal cardiac myocytes at physiological pacing rates – stable APs.
Stern (1992) showed that graded release is achieved when it is assumed that LCCs can only trigger Ca2+ release from their adjacent RyRs in the dyad. Under this assumption, graded release arises as the result of statistical recruitment of release clusters, a process known as local control of Ca2+ release. Guided by this insight, Greenstein and Winslow (2002) showed that when local control is incorporated into ventricular myocyte models by simulating the stochastic gating of LCCs and RyRs in each dyad, AP properties are stabilized. However, one drawback of models based on systems of stochastic ordinary differential equations is that solution of these equations is computationally demanding. Hinch et al. (2004) resolved this problem by using the fact that the time rate of change of dyadic Ca2+ concentration is so fast relative to the time evolution of any other biological process in the models that it can be assumed to immediately reach its steady-state value. This simple, reasonable assumption enabled the graded release property to be modeled using a low-dimensional system of ordinary differential equations in which LCCs and RyRs behave as a strongly coupled system. Incorporation of this “coupled LCC-RyR model” into cardiac ventricular myocyte models enabled these models to achieve graded release with high gain and stable APs (Greenstein et al., 2006) in a more computationally efficient manner. The advantage of such models, as compared to models with phenomenological formulations of the release mechanism, is that they can be used to study the functional consequences of altered molecular function on ECC gain since this property emerges as a result of capturing fundamental biological detail. This is not true of phenomenological models formulated using ECC gain functions that are explicitly built into the models.
In 2003, Cortassa et al. (2003) formulated a computational model of cardiac mitochondria including descriptions of the tri-carboxylic acid (TCA) cycle and its regulation by Ca2+, oxidative phosphorylation, the F1-F0 ATPase, the adenine nucleotide translocator, the Ca2+ uniporter, the Na+-Ca2+ exchanger, and mitochondrial Ca2+ dynamics. In 2006, this model was integrated into a version of the Jafri–Rice–Winslow model of the guinea pig ventricular myocyte (Jafri et al., 1998) that had been extended to include a description of isometric force generation (Rice et al., 2000). This integrative ECC/mitochondrial energetics (ECME) model (Cortassa et al., 2006) also described coupling between mitochondrial ATP production and energy requiring membrane transporters, as well as control of mitochondrial energetics by cytosolic Ca2+. This model was able to reconstruct steady-state relationships between force generation and oxygen consumption at different stimulus frequencies, as well as rapid temporal changes in mitochondrial NADH and Ca2+ in response to abrupt changes in workload. Nonetheless, this model is a common pool model exhibiting non-physiological all-or-none Ca2+ release. Incorporating the graded release property into this model is important because mitochondria are bounded at each end by the JSR Ca2+ release sites, a close association that is supported by the observation that there are electron dense structures linking the mitochondrial outer membrane to t-tubules (Hayashi et al., 2009). This possible colocalization of mitochondria and the Ca2+ release sites implies that mitochondria may sense the local dyadic Ca2+ signal rather than the bulk cytosolic Ca2+ signal exclusively. In addition, mitochondria are “buffers” of Ca2+ by virtue of the presence of the Ca2+ uniporter in the inner mitochondrial membrane. Therefore, mitochondria may not only sense and be regulated by the large, fast, local dyadic Ca2+ signal, they may also act to buffer this signal, thereby influencing ECC (Maack et al., 2006).
As a first step toward investigating these important questions, we present an extension to the ECME model incorporating the coupled LCC-RyR formulation of graded release and description of the local Ca2+ signal. We demonstrate that this model of the guinea pig ventricular myocyte is able to reconstruct a broad range of experimental data. The model predicts that interactions between voltage-dependent properties of ECC gain and AP shape during the plateau phase have an important role in the timing of the Ca2+ transient and thus force generation. This prediction, which emerges from the underlying graded release model, is validated by experimental data. Further, we show that factors influencing AP plateau shape such as magnitude of the fast transient outward K+ current (in species other than guinea pig) can significantly affect timing of Ca2+ release. This model prediction is also validated by experimental data. These behaviors are specific to the graded release model, and cannot be revealed when using common pool models with all-or-none release. Finally, we show preliminary results indicating that the model predicts experimentally measured effects of mitochondrial Ca2+ uniporter block on amplitude of the cytosolic and mitochondrial Ca2+ transients, demonstrating the important role of beat-to-beat Ca2+ buffering by mitochondria.
Materials and Methods
The Coupled LCC-RyR Ca2+ Release Unit
We have incorporated a coupled LCC-RyR model of CICR into the ECME guinea pig myocyte model of Cortassa et al. (2006; Figure 1A). The coupled LCC-RyR model of the Ca2+ release unit (CaRU) is based on that presented previously for canine myocytes by Greenstein et al. (2006; see Figures 1B,C). The CaRU is represented by a single LCC in the t-tubule membrane, a single RyR located in the closely opposed JSR membrane, and a dyadic volume in the space between them, which functions as a separate Ca2+ compartment (Figure 1A). The rate of Ca2+ diffusion from the dyadic space to the cytosol is sufficiently rapid allowing for the assumption that subspace Ca2+ levels equilibrate instantaneously and can therefore be expressed algebraically in terms of the fluxes through the LCC and RyR. Another simplification arises from the assumption that refilling of the JSR occurs quickly enough that the Ca2+ levels in the JSR can be assumed to be approximately equal to network SR (NSR) Ca2+ levels. In this minimal model, the single model RyR represents the estimated number of release channels per LCC measured in guinea pig (Bers and Stiffel, 1993), and thus corresponds to a cluster of five simultaneously gating RyRs. Therefore, unitary flux is increased to five times that of a single-channel. The CaRU model is made up of 40 states (Figure 1C), which represent all possible pairings of the 10 state LCC model and the four state RyR model. Further details on the coupled LCC-RyR formulation may be found in Greenstein et al. (2006).
Figure 2 demonstrates kinetic and steady-state properties of the LCC model. Parameters of the LCC model were constrained using voltage-clamp data obtained from isolated guinea pig ventricular myocytes measured at 34–37°C. Figure 2A shows the L-type Ca2+ current peak I–V relationship. ICa,L is non-zero for test potentials between approximately −40 and +60 mV with a maximal peak of −32 μA/μF at +10 mV. The membrane potential at which the I–V curve peaks is in good agreement with data from four different guinea pig studies (Rose et al., 1992; Allen, 1996; Grantham and Cannell, 1996; Linz and Meyer, 1998). Peak current is −32 μA/μF at +10 mV, at the high end of measured values. For comparison, experiments show −21 μA/μF at 37°C (Grantham and Cannell, 1996), −25.68 μA/μF temperature-adjusted from 34 (Allen, 1996) to 37°C, and −24 μA/μF temperature-adjusted from 22 (Rose et al., 1992; Allen, 1996; Grantham and Cannell, 1996; Linz and Meyer, 1998) to 37°C, where adjustments are made using a Q10 value of 2.96 from Cavalié et al. (1985) Steady-state CDI (Figure 2B), was constrained using data from double-pulse voltage-clamp protocols (Hadley and Lederer, 1991; Linz and Meyer, 1998), with CDI being greater than VDI at all potentials. VDI properties were constrained using data from Linz and Meyer (1998), who measured a non-specific current through LCCs in a Ca2+-free solution, and from Hadley and Lederer (1991), who determined VDI from measurements of LCC gating current charge immobilization. Rate of recovery from VDI was constrained using double-pulse voltage-clamp data from isolated rabbit ventricular myocytes (Mahajan et al., 2008). Figure 2C shows the time course of ICa,L at various test potentials. The current peaks 3 ms after stimulus before decaying over approximately 100 ms to a value near zero. The time course of ICa,L recordings, including time to peak ICa,L, are in good agreement with data from Linz and Meyer(1998; Figure 2D). Different peak magnitudes between model and experimental results suggest differing channel density, which is also reflected in Figure 2A. The number of LCCs (and thus release units) in the model was set to 339,000 in order to match experimental data on fractional release, as discussed in Section “CICR During the Action Potential” below. This number of LCCs is between the estimate of ≈276,000 predicted by binding experiments (Bers and Stiffel, 1993) and the estimate of ≈500,000 from LCC gating current studies (Hadley and Lederer, 1991).
Figure 2. Validation of the L-type Ca2+ current. (A) Current-voltage relation for the model (blue) compared with experimental data from Rose et al. (1992; green), Linz and Meyer(1998; red), Grantham and Cannell(1996; teal), and Allen (1996; purple). Recordings from Rose et al. (1992), Linz and Meyer (1998), and Allen (1996) were adjusted to 37°C using the Q10 value from Cavalié et al. (1985). (B) Steady-state availability in the presence of CDI and VDI (black) and with VDI only (blue) from the model (lines) is compared to experimental data from Linz and Meyer(1998; squares) and Hadley and Lederer(1991; triangles). (C) Model traces of ICa,L vs. time for 400 ms test pulses from a pre-clamp of −40 mV. (D) Experimental ICa,L traces from Linz and Meyer (1998), also from a pre-clamp of −40 mV. For the top set of traces in (C,D), the square represents a test potential of −10 mV and the triangle a test potential of 0 mV. For the bottom set, test potentials are 10, 30, and 50 mV for the circle, diamond, and triangle, respectively. (D) Was reproduced with copyright permission of John Wiley and Sons.
Ryanodine receptor properties are those from our previous model (Greenstein et al., 2006). Briefly, the RyR model comes from a formulation by Rice et al. (1999), which was modified from a model by Keizer and Smith(1998; Figure 1B). Upon elevation of subspace Ca2+ levels, RyRs rapidly transition from state 1 through state 2 into state 3, the open state. Termination of release occurs as the channel transitions from state 3 to state 4, the inactive state, where it remains until subspace Ca2+ levels drop and the RyR returns to state 1. Within the context of the whole-cell model, this RyR model produces an increasing load-dependent fractional release relationship. Figure 3 shows the fraction of total SR Ca2+ released by an AP at varying SR loads. The model relationship is qualitatively similar to experimental data from rabbit (Shannon et al., 2000) and ferret (Bassani et al., 1995), two species which have strong Na+-Ca2+ exchanger (NCX) contribution to relaxation, resulting in Ca2+-cycling that is more similar to guinea pig than rat or mouse.
Figure 3. Sarcoplasmic reticulum load-dependence of fractional release. Model SR load vs. fractional release (solid line) compared to experimental data from Shannon et al. (2000) in rabbit (dots) and Bassani et al. (1995) in ferret (triangles).
Ion Channels and Ca2+ Cycling
Models for the remaining (non-dyadic) channels, pumps, and exchangers are based on those of the 2006 ECME model (Cortassa et al., 2006) with the following changes. The previous formulation of the delayed rectifier potassium current IK is replaced with the formulation of Zeng et al. (1995) that separates this current into two components, the rapid and slow (IKr and IKs) delayed outward rectifier currents. IKr and IKs current amplitudes were adjusted to match model APD to experimental data (see Results). The equation for IKp was modified to better fit experimental data (Yue and Marban, 1988). The NCX model was updated with the model of Weber et al. (2001), which incorporates allosteric regulation by cytosolic Ca2+. NCX affinity constants for allosteric Ca2+ regulation and for intracellular Na+ were modified slightly to fit guinea pig NCX voltage dependant behavior (Maack et al., 2005) and diastolic Ca2+ levels under voltage-clamp (Han et al., 1994; Isenberg and Han, 1994). This allowed the model to maintain appropriate diastolic Ca2+ levels while still achieving a time constant of SR load rest decay of approximately 80 s (not shown), in qualitative agreement when compared with experimental data exhibiting a tau of 36 s (Terracciano et al., 1995). A mitochondrial Na+-H+ exchanger (Wei et al., 2011) was also added, which helps to ensure conservation of Na+ between the cytosolic and mitochondrial compartments. Ca2+ cycling was modified by adjusting NCX magnitude to achieve the approximately 60:40 balance of Ca2+ resequestration vs. Ca2+ export fluxes during diastole at a pacing frequency of 1 Hz (Bers, 2001). The ATP-dependent K+ current (IK(ATP)) was also added using a model by Ferrero et al. (1996) to incorporate the effects of ATP on APD. The conductance of IK(ATP) was modified to take into account the range of experimental data (Nichols and Lederer, 1990; Weiss et al., 1992). Minor changes were made to rate-limiting mitochondrial parameters in the TCA cycle, respiratory chain, and adenine nucleotide transporter to increase ATP supply at higher pacing frequencies. Ca2+-regulation of the mitochondria remains unchanged. Ca2+ uptake fluxes during 1 Hz pacing were numerically integrated at 1 ms resolution from the peak of the Ca2+ transient to the end of diastole (see Figure 4A) to simulate experimental relaxation protocols (Bers, 2001). Comparing the total moles of Ca2+ removed from the cytosol, the SR Ca2+-ATPase (SERCa) takes up 65.9% of the transported cytosolic Ca2+, NCX removes 28.9%, and the sarcolemmal (SL) Ca-pump removes 5.1%, close to estimates of 67% SERCa, 30% NCX, and 3% slow processes (Bers, 2001). The mitochondrial uniporter removes Ca2+ from the cytosol most effectively during the upstroke of the Ca2+ transient, while Ca2+ extrusion by the mitochondrial Na+/Ca2+ exchanger dominates during relaxation (see Figure 4B). The result is that a bolus of Ca2+ 12.7% of the size of the net integrated Ca2+ flux is directed back into the cytosol during relaxation. This illustrates the buffering role that mitochondria in this model play during Ca2+ cycling. Recirculation fraction, calculated as the amount of Ca2+ contributing to the Ca2+ transient that is derived from the SR stores at steady-state, is 69%. This is close to estimates between 63 and 67% (Bers, 2001). At steady-state, the amount of Ca2+ released from the SR must equal the amount resequestered during a given beat (see Figure 4B). Analogously the amount of Ca2+ entering the cell across the cell membrane must be extruded during relaxation. Thus the ratio of cytosolic Ca2+ extruded during relaxation to the amount resequestered into the SR is a good approximation of the ratio of Ca2+ flux into the cell across the sarcolemma to SR Ca2+ release. Here we have approximately 31% of cytosolic Ca2+ due to transsarcolemmal influx and 69% due to SR release. As expected at steady-state, total and mitochondrial integrated fluxes start and end at zero (see Figure 4B). Due to conservation of Ca2+, at steady-state the integrals of inbound and outbound fluxes across the sarcolemma also sum to zero.
Figure 4. Integrated cytosolic Ca2+ fluxes. Cytosolic Ca2+ uptake is given as the sum of SERCa, NCX, and sarcolemmal Ca2+-pump fluxes (A). SERCa contributes 65.9% of uptake, NCX 28.9%, and SL Ca-pump 5.1%. In this model mitochondria contribute significant beat-to-beat buffering. The mitochondrial uniporter predominates during the first 300 ms of the cycle, resulting in a net uptake of Ca2+ from the cytosol (B). For the remainder of the cycle, mitochondrial Na+-Ca2+ exchanger predominates, resulting in Ca2+ being transferred back out of the mitochondria to the cytosol. The mitochondria release an amount of Ca2+ equal to 12.7% of the net integrated Ca2+ uptake flux into the cytosol during relaxation. Integrating over the entire 1000 ms period (B), fluxes into the cytosol are given as positive and fluxes removing Ca2+ from the cytosol are negative. Fluxes shown represent SERCa (blue), NCX (red), SL Ca-pump (green), mitochondria (teal), background Ca2+ current (magenta), L-type Ca2+ current (yellow), SR release (gray), and total (black).
Model code was written in C++ using the SUNDIALS CVODE integration library and run using an IBM PC workstation with a 2.80-GHz processor and 2.5 GB of RAM. On such a workstation, this 75-state model takes approximately 15 s to compute 10 1 Hz APs with no output written to files. This is comparable to the performance of current models for which code is available (Shannon et al., 2004; Faber et al., 2007). Code and additional information are available online (http://www.icm.jhu.edu/models/). Equations and parameters for the model are also available in the online supplement.
CICR during the Action Potential
Figure 5A shows the model AP, which has a steady-state duration of 189 ms at 1 Hz pacing. This is consistent with experimental recordings obtained near the physiological temperature of 37°C with a 1-Hz APD range of approximately 130–180 ms (Arreola et al., 1991; Sicouri et al., 1996; Chen et al., 2000). AP shape, including resting membrane potential, peak voltage, and phase 2 and 3 AP slope, are in good correspondence with these experimental data.
Figure 5. Action potential and Ca2+ transient. Steady-state AP at 1 Hz pacing from the model (A) and experiment (Sipido et al., 1995b) (B) along with the corresponding Ca2+ transients (D,E). Note that the peak of the model transient is aligned with the middle of the AP plateau phase (dashed line), approximately 127 ms after stimulus. Experimental data in (E) show a similar delay in the Ca2+ transient peak of approximately 190 ms. L-type Ca2+ flux (JLCC) (C) shows a large, but brief peak aligned with the initial AP depolarization followed by a slow peak during the AP plateau. RyR flux (F) increases slowly and reaches its maximum in parallel with the JLCC slow peak and AP plateau. Dashed lines in (C,F) correspond to the time of Ca2+ transient peak from (D). Flux measurements are given with respect to subspace volume. (B,E) Were reproduced with copyright permission from The Physiological Society.
The model Ca2+ transient begins to increase as soon as ICa,L is triggered, but the peak is delayed 119 ms with respect to the peak voltage of the AP (Figures 5A,D). A delay of 144–190 ms is supported by experimental recordings (Sipido et al., 1995b; Grantham and Cannell, 1996; see also Figures 5B,E) and by the prediction of local control theory that local fidelity of ECC decreases at highly depolarized potentials, such as those that occur early in the guinea pig AP. At highly depolarized potentials, the driving force for Ca2+ current through LCCs is reduced, even though whole-cell LCC flux may be large (Figure 5C). The resulting smaller unitary currents are less likely to trigger release via RyRs, and hence ECC gain is reduced. This counteracts the increase in gain resulting from increased LCC open probability at more depolarized potentials. The model predicts that triggering of SR Ca2+ release is lowest during the early phases of the AP, and increases slowly as the AP evolves (Figure 5F) because of this initial low gain. This leads to a relative delay in the peak of the Ca2+ transient waveform as compared to other species. In contrast, the Greenstein et al. (2006) canine model, which also features the coupled LCC-RyR formulation of CICR, exhibits a cytosolic Ca2+ transient that peaks approximately 80 ms earlier. The difference can be attributed to the presence of an early repolarization phase in the canine AP, which quickly hyperpolarizes the membrane potential into the range where ECC gain is maximal. The guinea pig models of Faber et al. (2007) and Gaur and Rudy (2011), which represent the same species as the present model, but use different release formulations, also demonstrate earlier cytosolic Ca2+ peaks than shown in this model (see Discussion for a full analysis).
Figure 6 shows the characteristic fast peak, early decay, and late peak shape of the guinea pig model ICa,L, which is in good correspondence with experimental recordings (Arreola et al., 1991; Grantham and Cannell, 1996). Following the early decay, the current magnitude recovers partially in response to the increased driving force as the AP voltage drops, but the extent of this recovery is limited by increased Ca2+ release from the SR causing CDI. As the membrane potential repolarizes, the L-type current undergoes both VDI, as well as CDI resulting from SR Ca2+ release (see Figure 2B). LCC availability during the AP is shown in Figure 7. This quantity represents the occupancy of non-inactivated states in the channel model. The relation of LCC availability calculated from the model to membrane potential follows closely that of Linz and Meyer (1998). The kinetics vary between model and experiment due to the long APD in the experimental data. The APD90 from the experimental data is 380 ms, which more closely resembles the APD of a midmyocardial cell (Szigligeti et al., 1996), while the model is more representative of an epicardial cell. Simulations were also performed using the AP recording from the Linz and Meyer (1998) data as a voltage-clamp (Figures 7C,F). These simulations also show a rapid loss of availability, decreasing to a minimum of approximately 5% during the AP and a slow recovery of availability during the end of the repolarization phase.
Figure 6. ICa,L during the AP. Simulated ICa,L trace during steady-state AP at 1 Hz pacing (black) and experimental ICa,L trace from a guinea pig ventricular myocyte undergoing 1 Hz pacing at 35–37°C (blue). Blue trace is reproduced from Grantham and Cannell (1996) by copyright permission of the American Heart Association.
Figure 7. Inactivation of the L-type Ca2+ current during the AP. (A) Model 1 Hz AP. (D) Simulated LCC availability during the AP shown in (A). Experimental data from Linz and Meyer (1998) showing the AP clamp waveform (B) and LCC availability during the AP clamp (E) for guinea pig myocytes at 35°C. (C) Shows model membrane potential output using the trace from (B) as an AP clamp. The LCC availability during the AP clamp is given in (F). (B,E) Were reproduced by copyright permission of John Wiley and Sons.
The purpose of introducing the local control CaRU model is to incorporate a biophysically realistic representation of graded release in the myocyte model. As shown in Figure 8A, the RyR release and LCC trigger fluxes are smooth, continuous functions of membrane potential. The RyR release flux exhibits a peak that is displaced in the hyperpolarizing direction from that of the LCC trigger flux, a key feature of experimental measurements of Ca2+ release. This is again seen in the normalized RyR release flux (Figure 8B), better illustrating the representative leftward shift with relation to the LCC flux, as seen in the experimental data from Wier et al. (1994). Here, the LCC flux peaks at +10 mV while the RyR flux peaks at +5 mV. It is this shift that results in the characteristic monotonically decreasing gain function (Figure 8C). Voltage-dependent gain has not been measured in guinea pig, but the gain at 0 mV is close to experimental estimates in guinea pig (Sipido and Wier, 1991; Rocchetti et al., 2005).
Figure 8. Voltage-dependence of flux through LCCs and RyRs and ECC gain. (A) Voltage-dependence of maximal Ca2+ flux through LCCs (blue) and RyRs (green). (B) Normalized fluxes from (A). (C) ECC gain, as formulated by the ratio of maximal LCC flux to maximal RyR flux.
At 1 Hz pacing, the fractional release of total Ca2+ from the SR is approximately 33% (not shown). This quantity was calculated as unity minus the ratio of total systolic SR Ca2+ to total diastolic SR Ca2+. The 33% measurement given by the model agrees with 35% estimated in ferret (Bassani et al., 1995), a species with a similar recirculation fractional as the guinea pig (Bers, 2001). The peak RyR open probability is 2.9% at 1 Hz. As described by Bers (2001), given an experimentally measured peak release flux of 3 mM/s (Wier et al., 1994) and a unitary RyR flux near 0.4 pA (Mejía-Alvarez et al., 1999), only 40,000 RyRs need to be open at the time of maximal release. This is only 2% of the number of total RyRs calculated by Bers and Stiffel (1993).
Incorporation of the local control model of the CaRU into the myocyte model allows for the prediction of localized subspace Ca2+ levels (Figure 9). The model calculates subspace Ca2+ for four different dyad macrostates: with LCC and RyRs closed, with only the LCC open, with only the RyRs open, and with both the LCC and RyRs open. Average subspace Ca2+ can be estimated by summing over the predicted Ca2+ concentrations for these four scenarios, weighted by their respective probability of occurrence. During 1 Hz pacing, the predicted average subspace Ca2+ level peaks near 2 μM, four times higher than the peak of the cytosolic transient. Subspace Ca2+ for dyads with open LCCs and RyRs reaches a maximum of 45 μM during the AP plateau.
Figure 9. Predicted subspace Ca2+ levels. Model cytosolic Ca2+ transient during a steady-state 1 Hz AP (blue) and subspace Ca2+ transient averaged across all dyads (green). While the average subspace Ca2+ is approximately four times higher than that of the cytosol, the maximum subspace Ca2+ for a single dyad may reach 45 μM, measured as the maximum subspace Ca2+ for the open-open LCC-RyR configuration during a release event.
Action potential duration restitution describes the electrical response of the myocyte to a premature stimulus. When a myocyte is paced to steady-state at a constant basic cycle length (BCL), APD becomes constant from beat to beat. The time between the end of the AP and the onset of the next stimulus is the diastolic interval (DI). Immediately following an AP the cell will be in an inexcitable refractory state. As the DI increases beyond this refractory period toward the steady-state DI, an AP can be triggered, but its duration is less than that of the steady-state AP due to incomplete recovery from inactivation of ICa,L and INa. At the tissue level, slow APD restitution and high pacing frequency can lead to the formation of a functional conduction block as the wave of depolarization catches up with refractory tissue from the previous beat. This block fails to excite and changes the direction of wavefront propagation. Additionally, the relationship between DI and APD has been shown to influence whether such aberrant conduction patterns damp out or devolve into arrhythmias (Qu et al., 1999; Garfinkel et al., 2000).
An electrical restitution curve was generated for the model by pacing to steady-state at 2000 ms BCL and then saving the state of the model at the end of an AP. Using the state values from this time as initial conditions, premature stimuli at increasing DIs were applied and the resulting APD90 was measured. The results are compared to recordings under the same protocol from (midmyocardial) strips of guinea pig ventricle at 37°C (Figure 10A). A single exponential fit to the model results (yellow curve in Figure 10A) yields a time constant of 165 ms. This differs from measurements by Sicouri et al. (1996) and Davey et al. (2001) which are 26.7 and 78.4 ± 12.1 ms respectively for single guinea pig myocytes at 32°C. Possible causes of this discrepancy are presented in Section “Critique of the Model” below.
Figure 10. Restitution of APD. (A) APD restitution curves from guinea pig experimental data and models. Blue dots show 2000 ms BCL data from Bjornstad et al. (1993) connected by the biexponential fit reported by those authors. Green, red, and teal respectively show epicardial, midmyocardial, and endocardial experimental data from Sicouri et al. (1996). The purple line depicts output from the LRd07 model (Faber et al., 2007) and yellow is from the present model. A standard S1–S2 protocol was used to simulate APD restitution, starting from 0.5 Hz steady-state as in the protocol used by Sicouri et al. (1996). The single exponential time constant for restitution of the present model is 165 ms, compared with approximately 40 ms for Bjornstad et al. (1993) and approximately 46, 27, and 31 ms for Sicouri et al. (1996) epicardial, midmyocardial, and endocardial, respectively. (B) Alteration of the model restitution curve by reverting to the unmodified IKs formulation from Zeng et al. (1995). Blue is LRd07 model, green is the present model. Solid lines are simulations using the latest IKs formulation (one fast gate, one slow), dashed lines are for simulations with 1995 IKs formulation (two fast gates). Curves are normalized to the 2000-ms steady-state APD90 for each model, respectively.
Frequency-Dependence of APD and ECC
To determine the frequency-dependence of APD and accumulated force, the pacing protocol of Szigligeti et al. (1996) was followed. Briefly, the model was first paced to steady-state at 3000 ms BCL. Pacing frequency was then increased in a stepwise manner. APD, peak cytosolic Ca2+ and accumulated force were recorded after 3 min at 3000, 2000, 1500, 1000, 700, 500, and 300 ms BCL.
The model results presented in Figure 11A (blue line) show a decrease in APD with increasing pacing frequency, from an APD90 of 261 ms at 3000 ms BCL to an APD90 of 114 ms at 300 ms BCL. The rate of decrease of APD with increasing frequency follows closely the data of Szigligeti et al. (1996; Figure 11A; green line). Experimental studies have shown that intracellular sodium ([Na+]i) varies as a function of pacing frequency (Cohen et al., 1982; Wang et al., 1988; Maier et al., 1997). The model exhibits similar frequency-dependence of Na+ levels with average [Na+]i near 5.5 mM at 3000 ms BCL increasing to 8 mM by the end of the Szigligeti protocol. This Na+ accumulation is mediated through a frequency-dependent change in the balance of the Na+ entry via the fast sodium current (INa) that initiates the AP and the export of Na+ by the Na+-K+ pump during diastole. In this model, NCX plays a large role in early repolarization of the plateau at slow frequencies. As [Na+]i continues to decrease with an increase in BCL, outward NCX current is reduced, allowing prolongation of the AP plateau phase.
Figure 11. Frequency-dependence of APD and Force. (A) Comparison of model (blue) dependence of APD on frequency to that of Szigligeti et al. (1996; green). (B) Comparison of model (blue) dependence of force on frequency to that of Szigligeti et al. (1996; green). See Section “Frequency-Dependence of APD and ECC” of the text for conversion of normalized model output to force units.
The force model used here is the same as that implemented previously in the ECME model (Cortassa et al., 2006). The frequency-dependence of peak force follows closely the data from Szigligeti et al. (1996) in guinea pig papillary muscle (Figure 11B). As BCL decreases, isometric force increases exponentially as a result of the underlying increase in peak cytosolic Ca2+ concentration associated with increased SR load and gain. To compare the normalized force output of the model with those data the cross-sectional area used is 0.013 mm2, which is equivalent to a fiber with an elliptical profile of width 0.24 mm and thickness 0.069 mm, the minimum of the range of measurements for the Cortassa et al. (2006) right ventricular trabecula samples. The frequency-dependence of both cytosolic Ca2+ and isometric force play an important part in the ability of the model to reproduce experimental results on energy supply during work transitions (see below). Both cytosolic Ca2+ and force transients show faster peaks at higher pacing frequencies, as has been shown experimentally (Layland and Kentish, 1999; Bluhm et al., 2000; Davey et al., 2001; Raman et al., 2006; not shown). Cytosolic Ca2+ transients also show frequency-dependent acceleration of relaxation (not shown), such that the time of decay to half peak amplitude is decreased with increasing frequency. Force transients have been shown to have a similar property (Bers, 2001; Davey et al., 2001), but the rate of force transient relaxation in the model does not vary significantly across pacing frequencies.
Control of mitochondrial energy production is mediated via the Ca2+ sensitivity of key enzymes in the tri-carboxylic acid (TCA) cycle, and through regulation of the F1-F0 ATPase by ADP. As pacing frequency increases, a higher ADP:ATP ratio results from the increased ATP consumption at rapid contraction rates. Increased mitochondrial ADP levels stimulate the F1-F0 ATPase to generate ATP by utilizing the proton-motive force as an energy source. Flux of electrons through the electron transport chain would in and of itself deplete the NADH pool. However, the amplitude of the cytosolic Ca2+ transient also increases with pacing frequency, and this Ca2+ signal is communicated to the mitochondria via Ca2+ uptake by the mitochondrial Ca2+ uniporter. Elevated mitochondrial Ca2+ levels stimulate the TCA cycle to increase production of NADH to sustain a higher rate of respiration. As in the original ECME model (Cortassa et al., 2006), the rate of respiration for the model presented here increases with increasing pacing frequency. Respiration rate increases 3.3× between 4000 ms BCL pacing and 300 ms pacing, similar to Cortassa et al. (2006; their Figure 4A).
During abrupt changes in pacing frequency, the mechanism described above results in NADH transients with complex kinetics. While increased mitochondrial Ca2+ levels stimulate NADH supply, this occurs with slower kinetics than the increase in demand. The result is an abrupt decrease in NADH before restoration to a new steady-state at the higher pacing frequency. For decreases in pacing frequency, ATP demand drops causing an overshoot in NADH levels, since production is still stimulated by high mitochondrial Ca2+. Cytosolic Ca2+ transient amplitudes drop and mitochondrial Ca2+ levels eventually follow as Ca2+ is pumped out of the mitochondria by the mitochondrial Na+-Ca2+ exchanger. After tens or hundreds of seconds NADH levels reach a steady-state corresponding to the slower pacing frequency.
Figures 12A–C shows the response of the model to changes in workload. To simulate changes in heart rate, the model is paced at 0.25 Hz for 100 s, then pacing frequency is increased to a high workload frequency for 200 s before being allowed to recover at 0.25 Hz for another 200 s. This protocol is repeated for high workload frequencies of 0.5, 1, 1.5, and 2 Hz. Simulated NADH levels (Figure 12A) show an undershoot upon initiation of the high frequency stimulation, followed by a recovery to higher levels. Upon return to 0.25 Hz pacing, NADH levels exhibit an overshoot before beginning recovery to the steady-state 0.25 Hz level. In the model, NADH recovery to a steady-state is slow, taking approximately 600 s. The simulation results in Figure 12A were performed using the same experimental protocol. As a result, model NADH levels do not reach a steady-state within 200 s. While the kinetics and waveform of the NADH transients in the model are qualitatively similar to experimental data (Figure 12D), the lack of quantitative correspondence between the magnitudes of under- and over-shoots in our model, and the data of Brandes and Bers (1999), also seen in the original ECME model (Cortassa et al., 2006), may arise from the effects of ADP compartmentation. The smaller transient effects seen in the model indicate that the changes in workload in the cytosol are not inducing sufficiently large changes in the mitochondrial ADP pool. Averaged force (Figure 12B) shows an increase during high workload frequencies, up to 16.4 mN/mm2 at the end of 2 Hz pacing. These increases are indicative of larger cytosolic Ca2+ transients (not shown) and are followed by a recovery to rest levels at approximately 0.9 mN/mm2 during 0.25 Hz pacing. Mitochondrial Ca2+ levels (Figure 12C) show minimal beat-to-beat variation with significant increases in Ca2+ occurring gradually during periods of high pacing frequency. This property of mitochondrial Ca2+ loading is analogous to mitochondria acting as a low-pass filter of cytosolic Ca2+ changes. While the amplitude of cytosolic Ca2+ transients increases gradually after an increase in pacing frequency, Ca2+ levels always return to a diastolic value between beats. In contrast, mitochondrial Ca2+ accumulates almost monotonically during periods of high pacing frequency with only minor beat-to-beat variation.
Figure 12. Frequency-dependence of NADH levels. (A) Simulation of NADH concentration for a pacing protocol consisting of 100 s at 0.25 Hz [see labels above (A)], 200 s at higher pacing frequencies of 0.5, 1.0, 1.5, and 2.0 Hz then 200 s recovery at 0.25 Hz. Each 500 s protocol is started from 0.25 Hz steady-state initial conditions. (B) Moving average of model force output using a 4000-ms window. Force ranges from approximately 0.9 mN/mm2 at 0.25 Hz to 16.4 mN/mm2 at the end of the 2-Hz pacing period. (C) Model mitochondrial Ca2+ concentration. (D) Experimental data from Brandes and Bers (1999) follow a similar protocol. NADH signal is normalized to 0.25 Hz level. (D) Was reproduced from Brandes and Bers (1999) by copyright permission of the Biophysical Society.
To further illustrate the role of the Ca2+ uniporter in conveying cytosolic Ca2+ signals to the mitochondria, Figure 13 shows a comparison of model results with experimental data for uniporter block at 1 Hz pacing. After 75% of mitochondrial Ca2+ uniporters are blocked in the model, the cytosolic Ca2+ transient peak increases 51%, similar to the data shown from Maack et al. (2006). Conversely, uniporter block results in a decrease of the mitochondrial Ca2+ transient of about 60%. These effects also demonstrate the significance of beat-to-beat buffering of cytosolic Ca2+ by the mitochondria, since enough Ca2+ is taken up to reduce the amplitude of the cytosolic Ca2+ transient by approximately one-third. This also supports the results shown earlier (Figure 4) of Ca2+ fluxes during Ca2+ transient relaxation, where mitochondria exported an amount of Ca2+ into the cytosol 12.7% the size of the net cytosolic Ca2+ export flux. However, it is important to note that while the model mitochondria exchange large quantities of Ca2+ with the cytosol, strong buffering of Ca2+ in the mitochondrial matrix limits beat-to-beat variation of mitochondria Ca2+ transients to approximately 2% of peak mitochondrial Ca2+ at 1 Hz.
Figure 13. Effect of mitochondrial uniporter block. (A) Model results show the effect of 75% block of the mitochondrial uniporter (simulated by reducing the parameter Vmuni to 25% of its control value). Cytosolic Ca2+ transient magnitude is increased and (B) mitochondrial Ca2+ transient magnitude is decreased. Experimental data (Maack et al., 2006) for cytosolic Ca2+ levels (C) and mitochondrial Ca2+ levels (D) with addition of 10 nM Ru, a blocker of the mitochondrial uniporter.
In this model, the influence of the mitochondria on force production and vice versa is mediated largely through changes in the concentrations of cytosolic species. As described above, the Ca2+ buffering properties of the mitochondria affect the amplitude of the cytosolic Ca2+ transient. This in turn modulates the amplitude of the force transient. The myofibrils influence the mitochondria via their influence on ADP levels. Isometric contraction is linked to ADP:ATP levels via the acto-myosin ATPase, which is responsible for the bulk of ATP hydrolysis in the contracting myocyte. Increases in workload lead to a greater consumption of ATP, which results in a rise in ADP levels. This signal is conveyed from the cytosol to the mitochondria via creatine kinase acting on the diffusible creatine phosphate pool. It is the time delay imposed by this lengthy signaling cascade that leads to the disconnect between NADH supply and subsequent complex NADH transients.
The work presented here describes a mathematical model which represents the electrophysiology, Ca2+ cycling, isometric force development and mitochondrial energetics of the guinea pig cardiac myocyte. The novel feature of this model is the incorporation of a previously developed (Hinch et al., 2004; Greenstein et al., 2006) biophysically based model of local control of SR Ca2+ release. An advantage of this model is that it captures the key features of CICR (gradedness of Ca2+ release, voltage-dependent ECC gain) and does not require computationally expensive stochastic simulations of large numbers of individual ion channels (Greenstein and Winslow, 2002; Gaur and Rudy, 2011). We have previously shown that the inclusion of such a local control mechanism, which exhibits graded SR Ca2+ release, is crucial for the stability of the AP (Greenstein and Winslow, 2002), given that CDI of the LCC is much stronger than VDI (Hadley and Hume, 1987; Hadley and Lederer, 1991; Sipido et al., 1995a; Linz and Meyer, 1998;Peterson et al., 1999, 2000). Without a mechanistic description of this mechanism, common pool models are unstable because the strong negative feedback on ICa,L via CDI resulting from regenerative RyR Ca2+ release into the common pool essentially switches LCC trigger flux off prematurely. While other models (Viswanathan et al., 1999; Faber and Rudy, 2000; ten Tusscher et al., 2004; Mahajan et al., 2008) have been able to reproduce experimental features of ECC by incorporating phenomenological descriptions of graded release, use of phenomenological models always presents the risk that model predictions may be less reliable. Therefore the key question is: what do we gain from the formulation of a model that incorporates a biophysically based description of local control of Ca2+ release? This and other issues are discussed in the following sections.
Local Control Model Predicts Effects of AP Shape on Calcium-Release
The results of this study demonstrate an important functional relationship between early phase AP morphology and the kinetic properties of the cytosolic Ca2+ and force transients. The guinea pig ventricular myocyte model presented here, which includes an implementation of the new local control CaRU model (Greenstein et al., 2006), predicts that as a consequence of the shape of the guinea pig ventricular myocyte AP during the plateau phase, the Ca2+ transient peaks during the late phase of the AP (127 ms after the stimulus current; see Figure 14). This correlation of delayed Ca2+ transient peak with slowly repolarizing AP plateau voltage is seen in the experimental data (Arreola et al., 1991; Sah et al., 2003; Nishizawa et al., 2009). Consequently, the force transient is also delayed, having a peak that occurs after the AP is repolarized (Figure 14). This crucially important feature of the relative timing of the Ca2+ transient cannot be reconstructed using a common pool model, such as the prior version of the guinea pig ECME model (Cortassa et al., 2006). The delay in release timing exhibited by the present model is caused by a relatively low driving force for ICa,L early in the AP. While many LCCs open in the first few millisecond following the upstroke of the AP, unitary currents are small and the amount of Ca2+ entering any particular dyad via an open LCC is usually insufficient to trigger regenerative RyR opening. As the AP repolarizes, the driving force for ICa,L increases, LCC openings provide larger unitary currents, and hence are more effective at promoting the opening of RyRs, and release is more effectively triggered. This behavior, which was not built into the model, emerges as a result of incorporating the deterministic local control model which captures graded release and variable ECC gain. This model result emphasizes the role of the plateau potential in the nature of SR release triggering.
Figure 14. Action potential shape causes delayed Ca2+ transient and force. Normalized output from the guinea pig model is used to compare the kinetics of the AP, [Ca]i transient, and force transient. As seen above, the [Ca]i transient peaks near the end of the AP plateau. The force transient peak is further delayed and occurs after almost full repolarization of the cell.
As a consequence, it is important to note that differences in AP morphology (Figure 15A) can result in very different trigger L-Type Ca2+ currents (Figure 15B), and therefore Ca2+ transients with very different timing (Figure 15C). For example, the previously published canine myocyte model incorporating this same CaRU formulation (Greenstein et al., 2006) produced a cytosolic Ca2+ transient that peaked in the early phase of the AP (44 ms after the stimulus current; see Figures 15A,C). While the canine AP has a strong phase 1 repolarization due to the presence of the Ito current, guinea pig ventricular myocytes lack Ito (Findlay, 2003; Zhabyeyev et al., 2004) and their APs repolarize more slowly during the time period following the upstroke. When the guinea pig AP is at voltages near +40 mV, the canine AP has already repolarized to approximately +10 mV. As evident from the ECC gain curve (see Figure 8C), the canine AP will repolarize to a voltage at which gain is large very rapidly and trigger a relatively strong and synchronized release forming the early Ca2+ transient. In contrast, the slow repolarization of the guinea pig AP leads to a very slow increase in gain such that Ca2+ release is triggered in a more gradual manner, leading to a delay in the peak of the Ca2+ transient.
Figure 15. Impact of guinea pig and canine AP morphology on ICa,L and [Ca]i transients. In all panels blue traces are guinea pig model output and green traces are canine model output. (A) Comparison of APs from guinea pig and canine models (Greenstein et al., 2006). The canine AP has a significant early repolarization notch and a significantly longer APD. (B) L-type current traces peak near the same value, but guinea pig shows a much larger amount of late current. (C) Canine [Ca]i transient peak is approximately aligned with the AP notch, while the guinea pig [Ca]i transient peak occurs during the late plateau phase. (D) On the first beat after adding Ito,fast as in the Shannon et al. (2004) model with conductance of 0.2 mS/μF, the guinea pig AP exhibits a rapid initial repolarization and APD approaches that of canine. (E) With addition of Ito,fast, the guinea pig ICa,L trace exhibits a relation between fast peak and late current more similar to canine. The fast peak amplitude also increases substantially. (F) The peak of the guinea pig [Ca]i transient is now aligned with that of the canine model.
Many aspects of the canine and guinea pig local control models differ, including resting SR load, Ca2+ cycling parameters, and L-type Ca2+ channel parameters. To further test the hypothesis that AP shape is responsible for release timing, the control version of the guinea pig model was compared to a version with an AP which more closely resembles that of the canine model. The comparison version of the model incorporates a model of the Ito,fast current developed in the Shannon–Bers rabbit model (Shannon et al., 2004) representing the Ca2+-independent component of Ito (Ito,1) and consisting of one Hodgkin and Huxley (1952) activation gate and a second, slower Hodgkin–Huxley inactivation gate (Figures 15D–F). Both versions of the model were started from the control model’s 1 Hz initial conditions. The added Ito,fast model introduces a brief outward current upon depolarization to positive potentials. With this current included, the AP features rapid phase 1 repolarization and a Ca2+ transient that triggers with much less of a delay, similar to that of the canine model (Figures 15D–F). Addition of Ito,fast hyperpolarizes the initial membrane potential of the plateau to a voltage associated with high ECC gain. After LCCs open in response to the AP upstroke to approximately +50 mV, the Ito-induced repolarization reduces membrane potential and increases the driving force for ICa,L. This larger single-channel influx of Ca2+ triggers a larger Ca2+ flux from the RyRs by way of triggering release at a larger population of release sites than for the control AP at the same time. The result is that a larger release occurs earlier in the AP with Ito,fast, yielding a Ca2+ transient peak that occurs earlier in the AP than in control. Results were similar using an AP clamp recorded from the canine model, which yielded a Ca2+ transient peak of 1.05 μM occurring approximately 10 ms after the peak of the AP (not shown).
Use of a local control model such as this one featuring AP shape-dependent release will have important implications regarding behavior of tissue level model electro-mechanics. For example, in many species there are significant transmural differences in the expression levels of many key ion channels. The result is that APs from different tissue sites take on different morphologies (Nerbonne and Kass, 2005). This model predicts that these changes in AP shape will lead to changes in the timing of Ca2+ release, shape of the Ca2+ transient, and timing of force generation. AP clamp waveforms from guinea pig epicardial, midmyocardial, and endocardial cells as recorded by Sicouri et al. (1996) were input to the model to simulate the corresponding SR Ca2+ release and Ca2+ transient in each cell type (not shown). A 36% gradient of SERCa expression, measured with high resolution and accuracy by Anderson et al. (2011) was also applied. Epicardial cells were taken to have the control level of SERCa while SERCa expression was decreased by 18 and 36% in midmyocardial and endocardial cells, respectively. The resulting simulations show that release occurs earliest in midmyocardial cells due to the relatively hyperpolarized potential at the start of the AP plateau. Times to peak change by less than 10% without the application of a SERCa expression gradient, and rankings of times to peak remain the same between tissue types.
Given that this model predicts changes in release timing with AP morphology, the AP shape of the experimental species used may impact the strength of the conclusions that can be inferred about human electrophysiology. Among rabbit, canine, and human, all of which express Ito, the AP notch is more prominent in recordings from epicardial than endocardial myocytes (Fedida and Giles, 1991; Liu et al., 1993; Liu and Antzelevitch, 1995; Nabauer et al., 1996). For epicardial cells, canine APs exhibit a more pronounced notch than human, while rabbit APs have a noticeable early repolarization characteristic of Ito, but no strict spike-notch-dome morphology. Following phase 1 early repolarization, human AP plateau values peak below 25 mV for epicardial cells (Drouin et al., 1998; Piacentino et al., 2003). In comparison, canine plateau values for epicardial recordings at 1 Hz peak in the range of 10–15 mV (Liu et al., 1993; Liu and Antzelevitch, 1995) and rabbit plateau maxima take on values from 25 to 40 mV (Fedida and Giles, 1991; Puglisi et al., 1999). Guinea pig epicardial APs do not have a phase 1 repolarization and may take 50 ms for the plateau to decrease to 25 mV from the peak value near 50 mV (Sicouri et al., 1996). The current model predicts that these differences in notch depth and initial plateau height may significantly influence the timing of Ca2+ release and force generation in these different species. These results emphasize the importance of the inclusion of graded release in electromechanical models. Common pool electromechanical models without graded release are limited in their predictive scope. The all-or-none release produced by such common pool models fails to capture the sensitivity of the intracellular Ca2+ transient, and thus force transient, to changes in AP shape, as may occur with the reduction in Ito,1 in heart failure (Greenstein et al., 2006). Some of the earliest models incorporating electrophysiology with mechanical components were limited in their predictive scope because they modeled electrophysiology without any description SR release (Li et al., 2004) or with a description resulting in all-or-none release (Li et al., 2006). However as newer release descriptions succeeded in representing graded release (Chudin et al., 1999; Shannon et al., 2004; Mahajan et al., 2008), electromechanical models incorporating those advances gained the ability to predict mechanoelectrical feedback in arrhythmia (Jie et al., 2010; Keldermann et al., 2010) and patient-specific changes in electrical activation during heart failure (Aguado-Sierra et al., 2011).
In addition to AP morphology differences between species, diseases such as heart failure can produce significant alteration of AP morphology. In human heart failure Ito density has been shown to be downregulated (Beuckelmann et al., 1993; Nabauer et al., 1996). This is the major factor in the reduction of the AP notch in recordings from isolated failing cardiac myocytes. Another consequence to AP morphology is that the initial value of the plateau is elevated in heart failure. As discussed above, the membrane potential immediately following the AP peak plays a large role in the timing and magnitude of the Ca2+ transient. These findings corroborate those from Greenstein and Winslow (2002), which show that different Ito current densities have a significant impact on the efficacy of ECC, making AP morphology an important factor in predicting ECC properties during heart failure.
Critique of the Model
As with any computational model, compromises must be made in order to fit the range of experimental data for different protocols. The model restitution curves for 0.5 Hz pacing shown in Figure 10A gradually reach a plateau at which APD90 is approximately 219 ms. The time constant for this restitution (165 ms) is much slower than that observed in experiments (Bjornstad et al., 1993; Sicouri et al., 1996). This feature of the restitution curve in the model can be attributed to the time constant of IKs near the resting potential. Due to the slow deactivation time constants of IKs, there is a significant accumulation of active channels at DIs less than 1000 ms. The amount of current remains small because voltage is near the K+ reversal potential. If a premature stimulus shortens the DI, an even larger fraction of channels remain active, leading to larger IKs during the prematurely triggered AP. This results in a larger outward current that shortens the AP. Shorter DIs provide less time for IKs to deactivate, resulting in a larger IKs and shorter AP following the next stimulus. Thus the time constant of AP restitution is closely related to the time constant of IKs deactivation at negative potentials.
The model of IKs used here was formulated in the LRd99 model (Viswanathan et al., 1999) and was well fit to the two time constants of deactivation shown in experimental data (Matsuura et al., 1987). However, the Rudy group’s previous IKs model (Zeng et al., 1995), which shows a faster overall rate of deactivation, leads to faster restitution behavior that is a closer match to experimental data (Figures 10A,B), though still not as fast. Due to the dependence of IKs on the duration of diastole, this current’s kinetics affect both AP restitution and frequency-dependence behavior. The frequency-dependence of cytosolic Ca2+, force, and energy supply by mitochondria exhibited by the model is a good match to the experimental data (Szigligeti et al., 1996; Brandes and Bers, 1999). This model is unable to simultaneously achieve this frequency-dependent behavior and match the experimentally measured rate of AP restitution. Incorporation of a more detailed Markov model of IKs behavior (Silva and Rudy, 2005) may help reconcile these two cellular properties.
In addition, this model is not able to reproduce the Ca2+ restitution and related short-term interval-force relationships as described by Wier and Yue (1986). In other models (Rice et al., 2000; Faber et al., 2007) this behavior is mediated by a slow recovery of RyRs from inactivation. For short DIs, many RyRs remain inactivated so the maximal open probability is reduced, leading to smaller release and smaller Ca2+ transients. However, consistent with estimates by Bers (2001), this model has a much lower peak RyR open probability of 2.9% at 1 Hz pacing (compare to 100% for Faber et al., 2007 and 60% for Rice et al., 2000) and the proportion of channels that are closed and not inactivated does not drop below 70% at 1 Hz. The majority of RyRs are still available for activation in the event of a premature stimulus. As a result, the maximal RyR open probability triggered by a premature stimulus is not significantly reduced from control. This, coupled with the residual cytosolic Ca2+ stemming from the delayed Ca2+ transient exhibited by this model, limits the degree of Ca2+ restitution. A stimulus following a 1-Hz steady-state AP by a DI of 50 ms results in a premature AP with 26.5% shorter APD, but only a 7.7% smaller Ca2+ transient (not shown). In the coupled LCC-RyR formulation of Ca2+ release used in this model, JSR depletion cannot be tracked locally at individual release sites because the JSR and NSR are assumed to be in rapid equilibrium. Recent evidence suggests that the primary determinant of SR release termination and restitution of Ca2+ of sparks is the dynamics of local JSR refilling (Terentyev et al., 2002; Sobie et al., 2005; Zima et al., 2008b). Incorporation of a mechanistic formulation of this feature in future versions of this model should provide the ability to better reconstruct properties of Ca2+ restitution.
Computational models must weigh the advantages of physiological and mechanistic detail with computational efficiency in order to make useful predictions while still remaining tractable. Due to diffusion and compartmentation effects, the cardiac myocyte is subject to dynamic spatial gradients of a wide variety of ions and second messengers. In order to avoid the complexity of partial differential equations, the majority of models define compartments of uniform concentration. However the simple scenario of cytosolic, SR, mitochondrial, and dyadic compartments may not be sufficient to reproduce some experimental results. Experimental evidence (Weber et al., 2002, 2003) suggests that the concentrations of ions in close proximity to the sarcolemma may vary from those of the bulk cytosol. To account for this, several models feature a subsarcolemmal compartment (Shannon et al., 2004; Mahajan et al., 2008; Gaur and Rudy, 2011). While the model presented here was able to reproduce many critical ECC experiments without a submembrane volume, the addition of such a compartment may better account for the activity of sarcolemmal proteins that sense Na+ and Ca2+, such as NCX.
Greenstein and Winslow (2002) and Gaur and Rudy (2011) have previously published canine and guinea pig myocyte models, respectively, featuring a stochastic description of local control of SR Ca2+ release, a property that strongly motivated the work presented here. The ECC gain function of the low-dimensional ordinary differential equation model presented here is similar to that of the stochastic canine model. However, the ECC gain function of the Gaur and Rudy (2011) guinea pig model differs significantly from that shown here. Specifically, the gain in the model presented here decreases monotonically as membrane potential increases (Figure 8C) while the ECC gain of Gaur and Rudy (2011) decreases rapidly at hyperpolarized potentials and plateaus at a relatively high value at positive voltages. The gain for their model at +40 mV, which is the approximate peak of the AP, is roughly 10 times higher than that of the present model. As a result, Ca2+ release in their model (their Figure 3) is triggered fully at the start of the AP and terminates within 30 ms. In contrast, the lower gain in the present model at membrane potentials near the peak of the AP leads to incomplete, but persistent release that inactivates gradually during phase 3 repolarization. While, to the authors’ knowledge, voltage-dependent ECC gain has not been measured in guinea pig, the relative displacement of voltage-dependent LCC and RyR fluxes in this model closely resemble the shifted bell-shaped curves that give rise to the monotonic decreasing voltage-dependent ECC gain function measured in rat (Wier et al., 1994; Santana et al., 1996; Song et al., 2001).
An alternative approach to modeling graded release in deterministic myocyte models is to utilize more abstracted release descriptions. For example, the Rudy group (Faber et al., 2007; Gaur and Rudy, 2011) modeled the rate of Ca2+ flux through open RyRs as a function of the voltage-dependent Ca2+ trigger flux. Mahajan et al. (2008) utilized a simpler formulation lacking any description of RyR gating. Instead, the rate of Ca2+ release events was modeled as being proportional to the product of LCC open probability and a voltage-dependent gain function. One advantage of such models formulations is that they “break” the strong positive feedback loop that is well known to lead to all-or-none Ca2+ release in common pool models (Stern, 1992). These formulations can also reduce the number of model state variables, facilitating large-scale simulations at the tissue and whole-heart levels. In fact the Mahajan model (2008) is able to reproduce the late peak of the Ca2+ transient observed in this work through its phenomenological incorporation of voltage-dependent gain. However, a disadvantage of these more abstract models is that they cannot be used to predict the effects of events such as fundamental changes in RyR gating on ECC gain properties (Hashambhoy et al., 2010) without additional assumptions, because the gain function itself is built directly into the models.
The approach of approximating cellular function using more abstract (i.e., less mechanistic) descriptions to reduce the number of state variables could also be applied to the mitochondria. In a scenario where abnormal ATP supply does not need to be addressed, the Ca2+ buffering role of the mitochondria could be approximated by a slow buffer. Such an approximation would continue to influence the cytosolic Ca2+ transients on a beat-to-beat basis and over longer periods of stimulation by accumulating or releasing Ca2+ slowly as cytosolic Ca2+ peaks change. However, the slow buffer approximation would begin to break down under conditions where a large mitochondrial to cytosolic Ca2+ gradient is present, such as upon commencement of beta adrenergic stimulation.
The model presented here lacks a JSR compartment, which was combined with the network SR in developing the deterministic formulation of the CaRU model (Hinch et al., 2004). As a result, the model does not reproduce local SR depletion during release. Instead, this CaRU model relies on RyR inactivation for termination of release. While some experimental data supports inactivation as a possible mechanism of release termination (Sham et al., 1998), more recent work suggest that local JSR Ca2+ depletion may be the major mechanism for release termination (Terentyev et al., 2002; Zima et al., 2008a,b; Domeier et al., 2009; Stevens et al., 2009). The RyR model presented here precludes the possibility of incorporating local JSR depletion. In addition, other groups (Shannon et al., 2004; Mahajan et al., 2008; Gaur and Rudy, 2011) have incorporated SR load-dependence into RyR behavior, a feature not included in this model. Its inclusion may help reconstruct the steep SR load-dependent increase in fractional release seen in some datasets (Shannon et al., 2000).
Other models that incorporate both ECC and mitochondrial energetics exist. Magnus and Keizer (1998) developed one of the first comprehensive models of cytosolic and mitochondrial Ca2+ handling in pancreatic β-cells. Their Ca2+-regulation of the TCA cycle was the basis on which the energetics model used here was developed (Cortassa et al., 2003, 2006). Subsequently, Matsuoka et al. (2004) formulated a cardiac myocyte model describing electrophysiology and Ca2+ handling between mitochondria and cytosol. However, this model is not able to capture the biphasic response of NADH levels to changes in pacing frequency. Another recent model (Hatano et al., 2011) is novel in that it describes electrophysiology, mitochondrial energetics, and spatiotemporal changes in Ca2+ and many metabolites. While this model supports the existence of subsarcolemmal Ca2+ gradients as predicted by NCX calculations (Weber et al., 2002), and predicts the cellular distribution of ADP, the coarse resolution of the 3D mesh required for computational tractability precludes the ability to make predictions regarding subspace Ca2+ levels. The advantage of the model presented here is that its mechanistic description of ECC and Ca2+ cycling in the cell will provide a platform for future investigations in Ca2+-regulation of mitochondrial energy production.
We have developed a mechanistically detailed description of ECC in the guinea pig cardiac myocyte combined with modules describing energetics and isometric force. This model successfully reproduces key ECC properties of graded SR Ca2+ release and voltage-dependent gain. Additionally, the incorporation of mitochondrial energetics allows the model to qualitatively reproduce changes in NADH in response to changes in cardiac workload. Using this model we can improve our understanding of how changes in AP shape and Ca2+ transients affect energy supply and developed force in normal and failing myocytes.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This work was supported by the National Heart, Lung and Blood Institute (HL105216, HL105239) to R. L. Winslow and NSF and NDSEG graduate fellowships to L. D. Gauthier. R. L. Winslow is the Raj and Neera Singh Professor of Biomedical Engineering.
The Supplementary Material for this article can be found online at http://www.frontiersin.org/Computational_Physiology_and_Medicine/10.3389/fphys.2012.00244/abstract
Aguado-Sierra, J., Krishnamurthy, A., Villongco, C., Chuang, J., Howard, E., Gonzales, M. J., Omens, J., Krummen, D. E., Narayan, S., Kerckhoffs, R. C. P., and Mcculloch, A. D. (2011). Patient-specific modeling of dyssynchronous heart failure: a case study. Prog. Biophys. Mol. Biol. 107, 147–155.
Anderson, T., Wulfkuhle, J., Petricoin, E., and Winslow, R. L. (2011). High resolution mapping of the cardiac transmural proteome using reverse phase protein microarrays. Mol. Cell. Proteomics 10, M111.008037.
Arreola, J., Dirksen, R. T., Shieh, R. C., Williford, D. J., and Sheu, S. S. (1991). Ca2+ current and Ca2+ transients under action potential clamp in guinea pig ventricular myocytes. Am. J. Physiol. Cell Physiol. 261, C393–C397.
Bers, D. M., and Stiffel, V. M. (1993). Ratio of ryanodine to dihydropyridine receptors in cardiac and skeletal muscle and implications for E-C coupling. Am. J. Physiol. Cell Physiol. 264, C1587–C1593.
Bjornstad, H., Tande, P. L. M., Lathrop, D. A., and Refsum, H. (1993). Effects of temperature on cycle length dependent changes and restitution of action potential duration in guinea pig ventricular muscle. Cardiovasc. Res. 27, 946–950.
Bluhm, W. F., Kranias, E. G., Dillmann, W. H., and Meyer, M. (2000). Phospholamban: a major determinant of the cardiac force-frequency relationship. Am. J. Physiol. Heart Circ. Physiol. 278, H249–H255.
Cavalié, A., Mcdonald, T. F., Pelzer, D., and Trautwein, W. (1985). Temperature-induced transitory and steady-state changes in the calcium current of guinea pig ventricular myocytes. Pflugers Arch. 405, 294–296.
Chen, C. C., Lin, Y. C., Chen, S. A., Luk, H. N., Ding, P. Y. A., Chang, M. S., and Chiang, C. E. (2000). Shortening of cardiac action potentials in endotoxic shock in guinea pigs is caused by an increase in nitric oxide activity and activation of the adenosine triphosphate-sensitive potassium channel. Crit. Care Med. 28, 1713.
Cortassa, S., Aon, M. A., O’Rourke, B., Jacques, R., Tseng, H.-J., Marbán, E., and Winslow, R. L. (2006). A computational model integrating electrophysiology, contraction, and mitochondrial bioenergetics in the ventricular myocyte. Biophys. J. 91, 1564–1589.
Davey, P., Bryant, S., and Hart, G. (2001). Rate-dependent electrical, contractile and restitution properties of isolated left ventricular myocytes in guinea-pig hypertrophy. Acta Physiol. Scand. 171, 17–28.
Domeier, T. L., Blatter, L. A., and Zima, A. V. (2009). Alteration of sarcoplasmic reticulum Ca2+ release termination by ryanodine receptor sensitization and in heart failure. J. Physiol. (Lond.) 587, 5197–5209.
Faber, G. M., Silva, J., Livshitz, L., and Rudy, Y. (2007). Kinetic properties of the cardiac L-type Ca2+ channel and its role in myocyte electrophysiology: a theoretical investigation. Biophys. J. 92, 1522–1543.
Fabiato, A. (1985). Time and calcium dependence of activation and inactivation of calcium-induced release of calcium from the sarcoplasmic reticulum of a skinned canine cardiac Purkinje cell. J. Gen. Physiol. 85, 247–289.
Ferrero, J. M., Saiz, J., Ferrero, J. M., and Thakor, N. V. (1996). Simulation of action potentials from metabolically impaired cardiac myocytes: role of ATP-sensitive K+ current. Circ. Res. 79, 208–221.
Garfinkel, A., Kim, Y.-H., Voroshilovsky, O., Qu, Z., Kil, J. R., Lee, M.-H., Karagueuzian, H. S., Weiss, J. N., and Chen, P.-S. (2000). Preventing ventricular fibrillation by flattening cardiac restitution. Proc. Natl. Acad. Sci. U.S.A. 97, 6061–6066.
Hatano, A., Okada, J.-I., Washio, T., Hisada, T., and Sugiura, S. (2011). A three-dimensional simulation model of cardiomyocyte integrating excitation-contraction coupling and metabolism. Biophys. J. 101, 2601–2610.
Hayashi, T., Martone, M., Yu, Z., Thor, A., Doi, M., Holst, M., Ellisman, M., and Hoshijima, M. (2009). Three-dimensional electron microscopy reveals new details of membrane systems for Ca2+ signaling in the heart. J. Cell. Sci. 122, 1005–1013.
Hinch, R., Greenstein, J. L., Tanskanen, A. J., Xu, L., and Winslow, R. L. (2004). A simplified local control model of calcium-induced calcium release in cardiac ventricular myocytes. Biophys. J. 87, 3723–3736.
Liu, D.-W., and Antzelevitch, C. (1995). Characteristics of the delayed rectifier current (IKr and IKs) in canine ventricular epicardial, midmyocardial, and endocardial myocytes: a weaker IKs contributes to the longer action potential of the M cell. Circ. Res. 76, 351–365.
Liu, D. W., Gintant, G. A., and Antzelevitch, C. (1993). Ionic bases for electrophysiological distinctions among epicardial, midmyocardial, and endocardial myocytes from the free wall of the canine left ventricle. Circ. Res. 72, 671–687.
Maack, C., Ganesan, A., Sidor, A., and O’Rourke, B. (2005). Cardiac sodium-calcium exchanger is regulated by allosteric calcium and exchanger inhibitory peptide at distinct sites. Circ. Res. 96, 91–99.
Maack, C., Cortassa, S., Aon, M. A., Ganesan, A. N., Liu, T., and O’Rourke, B. (2006). Elevated cytosolic Na+ decreases mitochondrial Ca2+ uptake during excitation-contraction coupling and impairs energetic adaptation in cardiac myocytes. Circ. Res. 99, 172–182.
Mahajan, A., Shiferaw, Y., Sato, D., Baher, A., Olcese, R., Xie, L.-H., Yang, M.-J., Chen, P.-S., Restrepo, J. G., Karma, A., Garfinkel, A., Qu, Z., and Weiss, J. N. (2008). A rabbit ventricular action potential model replicating cardiac dynamics at rapid heart rates. Biophys. J. 94, 392–410.
Mejía-Alvarez, R., Kettlun, C., Ríos, E., Stern, M., and Fill, M. (1999). Unitary Ca2+ current through cardiac ryanodine receptor channels under quasi-physiological ionic conditions. J. Gen. Physiol. 113, 177–186.
Nabauer, M., Beuckelmann, D. J., Ãœberfuhr, P., and Steinbeck, G. (1996). Regional differences in current density and rate-dependent properties of the transient outward current in subepicardial and subendocardial myocytes of human left ventricle. Circulation 93, 168–177.
Nishizawa, H., Suzuki, T., Shioya, T., Nakazato, Y., Daida, H., and Kurebayashi, N. (2009). Causes of abnormal Ca2+ transients in Guinea pig pathophysiological ventricular muscle revealed by Ca2+ and action potential imaging at cellular level. PLoS ONE 4, e7069. doi:10.1371/journal.pone.0007069
Peterson, B. Z., Lee, J. S., Mulle, J. G., Wang, Y., De Leon, M., and Yue, D. T. (2000). Critical determinants of Ca2+-dependent inactivation within an EF-hand motif of L-type Ca2+ channels. Biophys. J. 78, 1906–1920.
Piacentino, V., Weber, C. R., Chen, X., Weisser-Thomas, J., Margulies, K. B., Bers, D. M., and Houser, S. R. (2003). Cellular basis of abnormal calcium transients of failing human ventricular myocytes. Circ. Res. 92, 651–658.
Puglisi, J. L., Yuan, W., Bassani, J. W. M., and Bers, D. M. (1999). Ca2+ influx through Ca2+ channels in rabbit ventricular myocytes during action potential clamp: influence of temperature. Circ. Res. 85, e7–e16.
Rocchetti, M., Besana, A., Mostacciuolo, G., Micheletti, R., Ferrari, P., Sarkozi, S., Szegedi, C., Jona, I., and Zaza, A. (2005). Modulation of sarcoplasmic reticulum function by Na+/K+ pump inhibitors with different toxicity: digoxin and PST2744 [(E,Z)-3-((2-Aminoethoxy)imino)androstane-6,17-dione hydrochloride]. J. Pharmacol. Exp. Ther. 313, 207–215.
Rose, W. C., Balke, C. W., Wier, W. G., and Marban, E. (1992). Macroscopic and unitary properties of physiological ion flux through L-type Ca2+ channels in guinea-pig heart cells. J. Physiol. (Lond.) 456, 267–284.
Sah, R., Ramirez, R. J., Oudit, G. Y., Gidrewicz, D., Trivieri, M. G., Zobel, C., and Backx, P. H. (2003). Regulation of cardiac excitation–contraction coupling by action potential repolarization: role of the transient outward potassium current (Ito). J. Physiol. (Lond.) 546, 5–18.
Santana, L. F., Cheng, H., Gómez, A. M., Cannell, M. B., and Lederer, W. J. (1996). Relation between the sarcolemmal Ca2+ current and Ca2+ sparks and local control theories for cardiac excitation-contraction coupling. Circ. Res. 78, 166–171.
Sham, J. S. K., Song, L.-S., Chen, Y., Deng, L.-H., Stern, M. D., Lakatta, E. G., and Cheng, H. (1998). Termination of Ca2+ release by a local inactivation of ryanodine receptors in cardiac myocytes. Proc. Natl. Acad. Sci. U.S.A. 95, 15096–15101.
Shannon, T. R., Ginsburg, K. S., and Bers, D. M. (2000). Potentiation of fractional sarcoplasmic reticulum calcium release by total and free intra-sarcoplasmic reticulum calcium concentration. Biophys. J. 78, 334–343.
Sipido, K. R., Callewaert, G., and Carmeliet, E. (1995a). Inhibition and rapid recovery of Ca2+ current during Ca2+ release from sarcoplasmic reticulum in Guinea pig ventricular myocytes. Circ. Res. 76, 102–109.
Sipido, K. R., Carmeliet, E., and Pappano, A. (1995b). Na+ current and Ca2+ release from the sarcoplasmic reticulum during action potentials in guinea-pig ventricular myocytes. J. Physiol. (Lond.) 489, 1.
Song, L.-S., Wang, S.-Q., Xiao, R.-P., Spurgeon, H., Lakatta, E. G., and Cheng, H. (2001). β-Adrenergic stimulation synchronizes intracellular Ca2+ release during excitation-contraction coupling in cardiac myocytes. Circ. Res. 88, 794–801.
Stevens, S. C. W., Terentyev, D., Kalyanasundaram, A., Periasamy, M., and Györke, S. (2009). Intra-sarcoplasmic reticulum Ca2+ oscillations are driven by dynamic regulation of ryanodine receptor function by luminal Ca2+ in cardiomyocytes. J. Physiol. (Lond.) 587, 4863–4872.
Szigligeti, P., Pankucsi, C., Bányász, T., Varró, A., and Nánási, P. P. (1996). Action potential duration and force-frequency relationship in isolated rabbit, guinea pig and rat cardiac muscle. J. Comp. Physiol. B Biochem. Syst. Environ. Physiol. 166, 150–155.
Terentyev, D., Viatchenko-Karpinski, S., Valdivia, H. H., Escobar, A. L., and Gyorke, S. (2002). Luminal Ca2+ controls termination and refractory behavior of Ca2+-induced Ca2+ release in cardiac myocytes. Circ. Res. 91, 414–420.
Terracciano, C. M. N., Naqvi, R. U., and Macleod, K. T. (1995). Effects of rest interval on the release of calcium from the sarcoplasmic reticulum in isolated Guinea pig ventricular myocytes. Circ. Res. 77, 354–360.
Weber, C. R., Ginsburg, K. S., Philipson, K. D., Shannon, T. R., and Bers, D. M. (2001). Allosteric regulation of Na/Ca exchange current by cytosolic Ca in intact cardiac myocytes. J. Gen. Physiol. 117, 119–132.
Wei, A.-C., Aon, M. A., O’Rourke, B., Winslow, R. L., and Cortassa, S. (2011). Mitochondrial energetics, pH regulation, and ion dynamics: a computational-experimental approach. Biophys. J. 100, 2894–2903.
Zeng, J., Laurita, K. R., Rosenbaum, D. S., and Rudy, Y. (1995). Two components of the delayed rectifier K+ current in ventricular myocytes of the Guinea pig type: theoretical formulation and their role in repolarization. Circ. Res. 77, 140–152.
Zhabyeyev, P., Asai, T., Missan, S., and Mcdonald, T. F. (2004). Transient outward current carried by inwardly rectifying K+ channels in guinea pig ventricular myocytes dialyzed with low-K+ solution. Am. J. Physiol. Cell Physiol. 287, C1396–C1403.
Zima, A. V., Picht, E., Bers, D. M., and Blatter, L. A. (2008a). Partial Inhibition of sarcoplasmic reticulum Ca release evokes long-lasting Ca release events in ventricular myocytes: role of luminal Ca in termination of Ca release. Biophys. J. 94, 1867–1879.
Keywords: calcium cycling, calcium-induced calcium-release, cardiac myocyte, computational model, excitation-contraction coupling, mitochondrial energetics
Citation: Gauthier LD, Greenstein JL and Winslow RL (2012) Toward an integrative computational model of the guinea pig cardiac myocyte. Front. Physio. 3:244. doi: 10.3389/fphys.2012.00244
Received: 15 March 2012; Paper pending published: 19 April 2012;
Accepted: 14 June 2012; Published online: 05 July 2012.
Edited by:Jennie Larkin, National Heart, Lung, and Blood Institute, USA
Reviewed by:Olga Solovyova, Institute Immunology and Physiology, Russia
Johan Elon Hake, Simula Research Laboratory, Norway
Copyright: © 2012 Gauthier, Greenstein and Winslow. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.
*Correspondence: Laura Doyle Gauthier, Department of Biomedical Engineering, Institute for Computational Medicine, Johns Hopkins University School of Medicine and Whiting School of Engineering, 3400 North Charles Street, Hackerman Hall 316A, Baltimore, MD 21218, USA. e-mail: firstname.lastname@example.org