The Role of Ca2+ Sparks in Force Frequency Relationships in Guinea Pig Ventricular Myocytes

Calcium sparks are the elementary Ca2+ release events in excitation-contraction coupling that underlie the Ca2+ transient. The frequency-dependent contractile force generated by cardiac myocytes depends upon the characteristics of the Ca2+ transients. A stochastic computational local control model of a guinea pig ventricular cardiomyocyte was developed, to gain insight into mechanisms of force-frequency relationship (FFR). This required the creation of a new three-state RyR2 model that reproduced the adaptive behavior of RyR2, in which the RyR2 channels transition into a different state when exposed to prolonged elevated subspace [Ca2+]. The model simulations agree with previous experimental and modeling studies on interval-force relations. Unlike previous common pool models, this local control model displayed stable action potential trains at 7 Hz. The duration and the amplitude of the [Ca2+]myo transients increase in pacing rates consistent with the experiments. The [Ca2+]myo transient reaches its peak value at 4 Hz and decreases afterward, consistent with experimental force-frequency curves. The model predicts, in agreement with previous modeling studies of Jafri and co-workers, diastolic sarcoplasmic reticulum, [Ca2+]sr, and RyR2 adaptation increase with the increased stimulation frequency, producing rising, rather than falling, amplitude of the myoplasmic [Ca2+] transients. However, the local control model also suggests that the reduction of the L-type Ca2+ current, with an increase in pacing frequency due to Ca2+-dependent inactivation, also plays a role in the negative slope of the FFR. In the simulations, the peak Ca2+ transient in the FFR correlated with the highest numbers of SR Ca2+ sparks: the larger average amplitudes of those sparks, and the longer duration of the Ca2+ sparks.


Introduction
Excitation-contraction (EC) coupling is the process triggered by the electrical excitation (depolarization) of the sarcolemma, through the Ca 2+ release from the sarcoplasmic reticulum (SR) to the contraction of the cardiac myocyte [1]. A frequency-dependent change in the cytoplasmic Ca 2+ transient underlies the force generated by the myocytes [2]. When external Ca 2+ enters the myocyte and causes the release of intracellular Ca 2+ from the SR, the elevation of Ca 2+ regulates the strength of the cardiac contraction. The change in contractile force of the cardiac myocytes at different pacing is termed an interval-force relationship or a force-frequency relationship (FFR). FFR is an essential intrinsic regulatory mechanism in cardiac myocyte contraction, to match the demand for increased blood supply [3,4].
The frequency-dependent contractile strength of the heart varies with the species of animal. In mammals, such as humans, rabbits, and guinea pigs, the relationship between cardiac contractile force and stimulation frequency (FFR) has been recorded to have a acting the L-type Ca 2+ channels in the dyadic subspace. The RyR2 model had four states: two closed states and two open states. Combining features of that model with the stochastic spark model [29] led to the development of a new three-state model: two closed states and one open state, as shown in Figure 1. The second closed state (C3) is an adaptive state. In the resting phase, almost all RyR2s stay in the closed state (C1); with the elevation of Ca 2+ in the dyadic subspace, the channels activate into an open state (O2) and, after some time, the channels might transition to an adapted state (C3). In this model, luminal regulation function (Ф), which depends on SR load ([Ca 2+ ]SR), modifies the channel opening rate. Therefore, the [Ca 2+ ]SR available to be released plays a major role in the developing force-frequency relationship [24]. The release of Ca 2+ from the SR is calculated by the equation: where is the maximal RyR2 single channel Ca 2+ Figure 2 shows the six-state L-type Ca 2+ channel model used in this work. In this model, states 2 (O2) and 3 (O3) are open, and states 1 (C1) and 6 (C6) are closed. C5 is the voltage-dependent inactivated state (VDI, O2  C5). C4, the Ca 2+ -dependent inactivated state (CDI, O2  C4), is controlled by the subspace Ca 2+ in each release site. Increased subspace Ca 2+ increases the rate of inactivation of LCC and prevents Ca 2+ overload in the myoplasm by limiting Ca 2+ entry [35]. The sixth state, C6, was added to the original fivestate model of Sun et al. [36], to allow for channels closure during the diastole. The CDI and VDI behaviors were constrained to match the experimental data from Morotti et al. [37]. In the resting phase, almost all RyR2s stay in the closed state (C 1 ); with the elevation of Ca 2+ in the dyadic subspace, the channels activate into an open state (O 2 ) and, after some time, the channels might transition to an adapted state (C 3 ). In this model, luminal regulation function (Φ), which depends on SR load ([Ca 2+ ] SR ), modifies the channel opening rate. Therefore, the [Ca 2+ ] SR available to be released plays a major role in the developing force-frequency relationship [24]. The release of Ca 2+ from the SR is calculated by the equation:

L-Type Ca 2+ Channel Model
where v 1 is the maximal RyR2 single channel Ca 2+ Figure 2 shows the six-state L-type Ca 2+ channel model used in this work. In this model, states 2 (O 2 ) and 3 (O 3 ) are open, and states 1 (C 1 ) and 6 (C 6 ) are closed. C 5 is the voltage-dependent inactivated state (VDI, O 2 → C 5 ). C 4 , the Ca 2+ -dependent inactivated state (CDI, O 2 → C 4 ), is controlled by the subspace Ca 2+ in each release site. Increased subspace Ca 2+ increases the rate of inactivation of LCC and prevents Ca 2+ overload in the myoplasm by limiting Ca 2+ entry [35]. The sixth state, C 6, was added to the original five-state model of Sun et al. [36], to allow for channels closure during the diastole. The CDI and VDI behaviors were constrained to match the experimental data from Morotti et al. [37].

L-Type Ca 2+ Channel Model
two closed states and two open states. Combining features of that model with the stochastic spark model [29] led to the development of a new three-state model: two closed states and one open state, as shown in Figure 1. The second closed state (C3) is an adaptive state. In the resting phase, almost all RyR2s stay in the closed state (C1); with the elevation of Ca 2+ in the dyadic subspace, the channels activate into an open state (O2) and, after some time, the channels might transition to an adapted state (C3). In this model, luminal regulation function (Ф), which depends on SR load ([Ca 2+ ]SR), modifies the channel opening rate. Therefore, the [Ca 2+ ]SR available to be released plays a major role in the developing force-frequency relationship [24]. The release of Ca 2+ from the SR is calculated by the equation:  Figure 2 shows the six-state L-type Ca 2+ channel model used in this work. In this model, states 2 (O2) and 3 (O3) are open, and states 1 (C1) and 6 (C6) are closed. C5 is the voltage-dependent inactivated state (VDI, O2  C5). C4, the Ca 2+ -dependent inactivated state (CDI, O2  C4), is controlled by the subspace Ca 2+ in each release site. Increased subspace Ca 2+ increases the rate of inactivation of LCC and prevents Ca 2+ overload in the myoplasm by limiting Ca 2+ entry [35]. The sixth state, C6, was added to the original fivestate model of Sun et al. [36], to allow for channels closure during the diastole. The CDI and VDI behaviors were constrained to match the experimental data from Morotti et al. [37]. During resting potential, all L-type channels are in a closed state (C 1 , and a change in the membrane potential activate them into an open state (O 2 ). A channel in O 2 state may continue to an open state (O 3 ) or transition into the voltage-dependent inactivated state (C 5 ), or excess Ca 2+ in dyadic subspace may bring them into the Ca 2+ -dependent inactivated state (C 4 ).

Simulation Frameworks
Simulations using the guinea pig ventricular myocyte model studied Ca 2+ transients in the cytosol at varying pacing frequencies. Force is assumed to be monotonically related to the size of the cytosolic calcium transient [38]. The foundation of the model was developed in 1 Hz pacing (basic cycle length (~1000 ms). All the parameters were adjusted for 1 Hz and compared all the plots with the experimental results and the previous models work. After this, the model simulated FFR, with basic cycle length as low as 0.2 Hz (1 beat in 5 s) to as high as 8 Hz (8 beats in 1 s). The other simulation frequencies are 0.20 Hz, 0.25 Hz, 0.33 Hz, 0.5 Hz, 1 Hz, 2 Hz, 3 Hz, 4 Hz, 5 Hz, 6 Hz, and 7 Hz. For every pacing frequency, each simulation ran for 10 s to reach a stable train of transients. The simulation data were recorded and plotted against the time. The amplitude and duration of AP, amplitudes of L-type current (I LCC ), Na + -Ca 2+ current (I ncx ), Na + current (I Na ), transient slow outward K + current (I ktos ), transient outward K + current (I ktof ), inward rectifier (I K1 ) K + current, cytoplasmic Ca 2+ concentration ([Ca 2+ ] myo ), NSR Ca 2+ concentration ([Ca 2+ ] nsr ), and RyR2 opening probability (P open ) from each beat were measured. Finally, the peak values of each ionic current were collected and plotted against the pacing frequencies.
The subcellular and molecular analysis of the FFR at Ca 2+ level was performed by calculating average numbers of Ca 2+ sparks, Ca 2+ spark amplitudes, and the Ca 2+ spark duration. Ten beats were used to compute these data after the simulation arrived in a stable state. After that, an analysis algorithm chose the segments of the systolic phase only and excluded any data from the diastolic phase. From the systolic segment, the dyadic Ca 2+ transient, which was at least 25 µM tall, was counted as a spark and used to sum the total sparks in that beat. Similarly, the sparks were counted from the rest of the beats, and the average numbers of the sparks were calculated. Likewise, the average amplitudes of the sparks from different beats of the same frequency were measured by adding the amplitude of each beat and dividing the total amplitude by the number of beats. After collecting average spark amplitudes of each beat, the mean amplitude was computed over 10 beats. The average spark duration was calculated using a similar protocol to that used for Ca 2+ amplitudes. The total spark durations were combined in each beat, the sum was divided by the number of sparks, and the mean spark duration was concluded from the 10 beats.

Numerical Methods
PGI CUDA Fortran was used to compile, execute, and simulate the program on a workstation using the Ubuntu Linux operating system. CUDA (compute unified device architecture) is a parallel computing programming language developed by NVIDIA for graphic processing units (GPUs). The workstations used for the simulations contain Fermibased C2050 graphics processing cards with CUDA Toolkit 6.0 and higher. To capture calcium dynamics at a single-channel level, a novel computational algorithm Ultra-Fast Markov chain Monte Carlo (UMCMC) method was used for the stochastic gating from the Ca 2+ release units [39]. Numerical integration of all ordinary differential equations used the explicit Euler method, with an adaptive time step that ranged from 10 ns to 100 ns.

Ca 2+ Transient FFR Curves
In the simulations, the amplitude of the Ca 2+ peak transient ([Ca 2+ ] myo ) rose with the increase in the pacing frequency, from 0.20 Hz (0.33 µM) to 4 Hz (0.64 µM), peaking at 4 Hz ( Figure 3A). The peak amplitude declines from 5 Hz (0.61 µM) to 7 Hz (0.56 µM). The model showed a positive slope until 4 Hz frequency, and a negative slope after that. Overall, as shown in Figure 3A, the force-frequency relation (FFR) curve was similar to experiments in guinea pig cardiac papillary muscles, as presented in Figure 3B [7,40].

Ca 2+ Transient FFR Curves
In the simulations, the amplitude of the Ca 2+ peak transient ([Ca 2+ ]myo) rose with the increase in the pacing frequency, from 0.20 Hz (0.33 µM) to 4 Hz (0.64 µM), peaking at 4 Hz ( Figure 3A). The peak amplitude declines from 5 Hz (0.61 µM) to 7 Hz (0.56 µM). The model showed a positive slope until 4 Hz frequency, and a negative slope after that. Overall, as shown in Figure 3A, the force-frequency relation (FFR) curve was similar to experiments in guinea pig cardiac papillary muscles, as presented in Figure 3B [7,40].  Figure 4D) behaved similarly to FFR, increasing up to 4 Hz pacing and decreasing thereafter. RyR2 adaptation, on the other hand, increased until 6 Hz and then declined ( Figure 4E). When the pacing frequency increases in the model, the APD decreased from 0.2 Hz to 7 Hz ( Figure 4F).  Figure 4D) behaved similarly to FFR, increasing up to 4 Hz pacing and decreasing thereafter. RyR2 adaptation, on the other hand, increased until 6 Hz and then declined ( Figure 4E). When the pacing frequency increases in the model, the APD decreased from 0.2 Hz to 7 Hz ( Figure 4F).

L-Type Current Decreases and I ncx Current Increases with the Rapid Pacing
The SR Ca 2+ concentration rises with the increasing pacing rate, while the Ca 2+ current in each beat via L-type channel decreases ( Figure 5A). On the other hand, the electrogenic I ncx current becomes elevated with the increase in the pacing frequency ( Figure 5B). Both these trends are due to the increase in the time-averaged myoplasmic [Ca 2+ ]. The Ca 2+dependent inactivation of L-type channels is attenuated at lower frequencies ( Figure 5C) and rises with the increase in frequency ( Figure 5D) (1 Hz vs. 6 Hz). The frequencydependent increase in RyR2 P open caused an increase in the SR Ca 2+ release, which in turn caused the LCC channels to become inactive ( Figure 5C,D). The loss of L-type current leads to a decrease in the plateau phase in AP, which also decreases the APD with an increasing pacing rate ( Figure 4F). The SR Ca 2+ serves as the feedback mechanism to the L-type channels, and their amplitude decreases with the increase in pacing frequency [41]. In this simulation, the extrusion of Ca 2+ from Na + -Ca 2+ exchanger constantly increases from 0.2 Hz to 7 Hz ( Figure 5B).

Calcium Sparks Are the Subcellular Mechanisms of FFR
The above results demonstrate the role of Ca 2+ transient in reproducing force-frequency relationship (FFR). The model allows an analysis of how Ca 2+ spark frequency and Ca 2+ spark amplitude regulate the Ca 2+ transient and the resulting contractile force of a myocyte. Ca 2+ spark properties show similar behavior to Ca 2+ transients ( Figure 3A). The average number of Ca 2+ sparks in each beat at different frequencies increases from 0.2 Hz to 4 Hz pacing, and the number of sparks gradually decreased thereafter ( Figure 6A). The maximum number of sparks appeared at 4 Hz pacing (83,553 ± 5105). The model displayed spark numbers which, when compared to 4 Hz value, decrease by~5%,~11%, 17% with 5, 6, and 7 Hz, respectively. Lukyanenko et al. [42] found that the frequency of the sparks increases with the increase in SR Ca 2+ load. The Ca 2+ spark frequency was very high at 4 Hz frequency. In addition to the spark frequency, the average spark amplitudes followed the same trend, with the peak amplitude (62.95 ± 0.55 µM) also occurring at 4 Hz pacing and decreasing at higher pacing frequencies. The SR load also affects the Ca 2+ spark amplitudes [43]. The model found the largest average amplitudes when SR load was higher at 4 Hz ( Figure 6B). A larger amplitude means a greater force would be generated at 4 Hz pacing. The Ca 2+ spark duration is another activity that showed the same kind of behavior, demonstrating the highest value at 4 Hz ( Figure 6C). The model has shown that the total Ca 2+ release was maximum at 4 Hz ( Figure 6D). The total Ca 2+ release was approximated by multiplying Ca 2+ spark amplitude and duration and number of sparks and dividing by 2, based on the assumption of a triangular Ca 2+ flux profile.

L-Type Current Decreases and Incx Current Increases with the Rapid Pacing
The SR Ca 2+ concentration rises with the increasing pacing rate, while the Ca 2+ current in each beat via L-type channel decreases ( Figure 5A). On the other hand, the electrogenic Incx current becomes elevated with the increase in the pacing frequency ( Figure  5B). Both these trends are due to the increase in the time-averaged myoplasmic [Ca 2+ ]. The Ca 2+ -dependent inactivation of L-type channels is attenuated at lower frequencies ( Figure 5C) and rises with the increase in frequency ( Figure 5D) (1 Hz vs. 6 Hz). The frequency-dependent increase in RyR2 Popen caused an increase in the SR Ca 2+ release, which in turn caused the LCC channels to become inactive ( Figure 5C,D). The loss of Ltype current leads to a decrease in the plateau phase in AP, which also decreases the APD with an increasing pacing rate ( Figure 4F). The SR Ca 2+ serves as the feedback mechanism to the L-type channels, and their amplitude decreases with the increase in pacing frequency [41]. In this simulation, the extrusion of Ca 2+ from Na + -Ca 2+ exchanger constantly increases from 0.2 Hz to 7 Hz ( Figure 5B).

Calcium Sparks Are the Subcellular Mechanisms of FFR
The above results demonstrate the role of Ca 2+ transient in reproducing for frequency relationship (FFR). The model allows an analysis of how Ca 2+ spark frequen and Ca 2+ spark amplitude regulate the Ca 2+ transient and the resulting contractile force a myocyte. Ca 2+ spark properties show similar behavior to Ca 2+ transients ( Figure 3A The average number of Ca 2+ sparks in each beat at different frequencies increases fro 0.2 Hz to 4 Hz pacing, and the number of sparks gradually decreased thereafter (Figu 6A). The maximum number of sparks appeared at 4 Hz pacing (83,553 ± 5105). The mo el displayed spark numbers which, when compared to 4 Hz value, decrease by ~5 ~11%, ~17% with 5, 6, and 7 Hz, respectively. Lukyanenko et al. [42] found that the f quency of the sparks increases with the increase in SR Ca 2+ load. The Ca 2+ spark freque cy was very high at 4 Hz frequency. In addition to the spark frequency, the avera spark amplitudes followed the same trend, with the peak amplitude (62.95 ± 0.55 µ also occurring at 4 Hz pacing and decreasing at higher pacing frequencies. The SR lo also affects the Ca 2+ spark amplitudes [43]. The model found the largest average amp tudes when SR load was higher at 4 Hz ( Figure 6B). A larger amplitude means a grea

Luminal Dependence and SR Ca 2+ Play Major Role in FFR
Luminal Ca 2+ regulates the RyR2 P open and thereby regulates the dynamics of intracellular Ca 2+ . Higher luminal Ca 2+ availability increases RyR2 open probability, while a lower value reduces it. Simulations showed that lowering luminal dependency by 20%, in comparison to the control luminal value, reduced the Ca 2+ transient peaks by~6% ( Figure 7A). The decrease in the luminal dependency lowers the RyR2 P open ( Figure 7B) and decreases the amount of Ca 2+ released from the SR Ca 2+ via the RyR2, thereby reducing CICR. The consequence of that drop is an increase in the SR Ca 2+ load, which is shown in Figure 7C. The decreasing luminal activity also played a role in lowering the diastolic fraction of RyR2, because a decrease in the cytosolic Ca 2+ lowered the RyR2 adaptation rate ( Figure 7D), which solely depends upon cytosolic Ca 2+ . Similarly, the numbers of Ca 2+ sparks ( Figure 7E) and the average spark amplitude ( Figure 7F) were lower with the reduced luminal Ca 2+, because the decreases in SR Ca 2+ load lowers the RyR2 P open and therefore decreases the number of Ca 2+ sparks. Simulations showed that 8% (57,837 ± 18,935 vs. 53,322 ± 17,377) fewer sparks per beat were released with a 20% reduction in luminal dependence than in the control simulation.
force would be generated at 4 Hz pacing. The Ca 2+ spark duration is another activity that showed the same kind of behavior, demonstrating the highest value at 4 Hz ( Figure 6C). The model has shown that the total Ca 2+ release was maximum at 4 Hz ( Figure 6D). The total Ca 2+ release was approximated by multiplying Ca 2+ spark amplitude and duration and number of sparks and dividing by 2, based on the assumption of a triangular Ca 2+ flux profile.

Luminal Dependence and SR Ca 2+ Play Major Role in FFR
Luminal Ca 2+ regulates the RyR2 Popen and thereby regulates the dynamics of intracellular Ca 2+ . Higher luminal Ca 2+ availability increases RyR2 open probability, while a lower value reduces it. Simulations showed that lowering luminal dependency by 20%, in comparison to the control luminal value, reduced the Ca 2+ transient peaks by ~6% ( Figure 7A). The decrease in the luminal dependency lowers the RyR2 Popen ( Figure 7B) and decreases the amount of Ca 2+ released from the SR Ca 2+ via the RyR2, thereby reducing CICR. The consequence of that drop is an increase in the SR Ca 2+ load, which is shown in Figure 7C. The decreasing luminal activity also played a role in lowering the diastolic fraction of RyR2, because a decrease in the cytosolic Ca 2+ lowered the RyR2 adaptation rate ( Figure 7D), which solely depends upon cytosolic Ca 2+ . Similarly, the numbers of Ca 2+ sparks ( Figure 7E) and the average spark amplitude ( Figure 7F) were lower with the reduced luminal Ca 2+, because the decreases in SR Ca 2+ load lowers the RyR2 Popen and therefore decreases the number of Ca 2+ sparks. Simulations showed that ~8% (57,837 ± 18,935 vs. 53,322 ± 17,377) fewer sparks per beat were released with a 20% reduction in luminal dependence than in the control simulation.

Adaptation Brings Negative Feedback Mechanism to the Open Probability of RyR2
In the model, adaptation shifts the modal gating behavior of RyR2. A 20% reduction of the adaptation rate of the RyR2 channel is accompanied by a 4% increase in the size of Ca 2+ transients ( Figure 8A), because the RyR2 P open ( Figure 8B) went up with the reduction of the RyR2 adaptation. The increased RyR2 P open increase SR Ca 2+ leak to lower peak diastolic SR Ca 2+ ( Figure 8C) while producing larger Ca 2+ transients. The RyR2 adaptation fraction decreased as expected ( Figure 8D). The increased SR release also produced a greater number of Ca 2+ sparks ( Figure 8E) and larger average spark amplitudes ( Figure 8F). In analyzing Ca 2+ spark behavior, a 5% increase (57,837 ± 18,937 vs. 60,768 ± 19,752) in the Ca 2+ sparks per beat was observed with this simulation.

The Role of RyR2 Opening Rate Constant in FFR
The opening probability plays a major role in the initiation of CICR. A 20% reduction in the RyR2 opening rate constant decreased the peak amplitude of Ca 2+ transient by 4% ( Figure 9A) compared with the control case. The decreased opening rate of RyR2 decreased RyR2 Popen ( Figure 9B), reducing SR Ca 2+ leak and thereby increasing diastolic SR load ( Figure 9C). The diastolic adaptation of RyR2 ( Figure 9D) remains almost the same, because the Ca 2+ transient did not increase sufficiently when compared to control condi-

The Role of RyR2 Opening Rate Constant in FFR
The opening probability plays a major role in the initiation of CICR. A 20% reduction in the RyR2 opening rate constant decreased the peak amplitude of Ca 2+ transient by 4% ( Figure 9A) compared with the control case. The decreased opening rate of RyR2 decreased RyR2 P open (Figure 9B), reducing SR Ca 2+ leak and thereby increasing diastolic SR load ( Figure 9C). The diastolic adaptation of RyR2 ( Figure 9D) remains almost the same, because the Ca 2+ transient did not increase sufficiently when compared to control conditions. The Ca 2+ sparks ( Figure 9E) and spark amplitude ( Figure 9F)

The Role of L-Type Ca 2+ Channels and NCX in FFR
The previous simulations showed that change properties of the RyR2 could affect FFR. The Ca 2+ flux through RyR2 play a critical role in SR Ca 2+ levels, as they regulate the efflux of Ca 2+ from the SR. On the other hand, the total cellular Ca 2+ is governed by the transport of Ca 2+ across the sarcolemma. When the NCX activity is increased by 25%, the FFR has a negative slope ( Figure 10A). This occurs because the increase Ca 2+ efflux across the sarcolemma limits the increase in SR Ca 2+ with increasing pacing frequency ( Figure 10B). In contrast, when Ca 2+ entry via the L-type Ca 2+ channels are reduced by 20%, negative FFR is observed in the simulations ( Figure 10C). The reduction in Ca 2+ influx reduced the filling of the SR with increasing pacing frequency, causing the negative FFR ( Figure 10D).

The Role of L-Type Ca 2+ Channels and NCX in FFR
The previous simulations showed that change properties of the RyR2 could affect FFR. The Ca 2+ flux through RyR2 play a critical role in SR Ca 2+ levels, as they regulate the efflux of Ca 2+ from the SR. On the other hand, the total cellular Ca 2+ is governed by the transport of Ca 2+ across the sarcolemma. When the NCX activity is increased by 25%, the FFR has a negative slope ( Figure 10A). This occurs because the increase Ca 2+ efflux across the sarcolemma limits the increase in SR Ca 2+ with increasing pacing frequency ( Figure 10B). In contrast, when Ca 2+ entry via the L-type Ca 2+ channels are reduced by 20%, negative FFR is observed in the simulations ( Figure 10C). The reduction in Ca 2+ influx reduced the filling of the SR with increasing pacing frequency, causing the negative FFR ( Figure 10D).

Pacing Protocols in FFR
Simulations changing the stimulation rate from 0.5 Hz pacing to 1.5 Hz pacing and back to 0.5 Hz pacing are shown in Figure 10. These simulations showed a minimal decrease in the amplitude of AP (39.49 to 39.23 mV) ( Figure 11A) from lower to higher and a 9 ms difference in the AP duration (172 ms to 163 ms). Similar to the classic staircase response, the myoplasmic Ca 2+ transient amplitude increases with a switch to 1.5 Hz pacing rate from 0.5 Hz pacing rate ( Figure 11B) and reaches a new steady state due to increased diastolic SR content ( Figure 11C). The fraction of RyR2 channels in the adapted state increases with the faster pacing rate, due to the higher Ca 2+ levels and shorter time for recovery before the subsequent action potential ( Figure 11D). On the other hand, when the frequency changes from a 1.5 Hz to 0.5 Hz, there is a transient increase in the Ca 2+ transient amplitude, followed by a decline to the steady state amplitude. This is due to the recovery for the RyR2 adaptation while the SR load is still elevated during the first set of beats after the rate change.

Pacing Protocols in FFR
Simulations changing the stimulation rate from 0.5 Hz pacing to 1.5 Hz pacing and back to 0.5 Hz pacing are shown in Figure 10. These simulations showed a minimal decrease in the amplitude of AP (39.49 to 39.23 mV) ( Figure 11A) from lower to higher and a 9 ms difference in the AP duration (172 ms to 163 ms). Similar to the classic staircase response, the myoplasmic Ca 2+ transient amplitude increases with a switch to 1.5 Hz pacing rate from 0.5 Hz pacing rate ( Figure 11B) and reaches a new steady state due to increased diastolic SR content ( Figure 11C). The fraction of RyR2 channels in the adapted state increases with the faster pacing rate, due to the higher Ca 2+ levels and shorter time for recovery before the subsequent action potential ( Figure 11D). On the other hand, when the frequency changes from a 1.5 Hz to 0.5 Hz, there is a transient increase in the Ca 2+ transient amplitude, followed by a decline to the steady state amplitude. This is due to the recovery for the RyR2 adaptation while the SR load is still elevated during the first set of beats after the rate change.

Discussion
The frequency-dependent performance of the myocardium changes the contractile strength of the heart. For humans, rabbits or guinea pigs, this force-frequency curve shows biphasic behavior, initially increasing and then decreasing [3,[5][6][7][8]. The force generated by the cardiac myocyte is a monotonically increasing function of the amplitude of the myoplasmic Ca 2+ transient [44]. A newly developed model of the ventricular myocytes of a guinea pig was used to study the role Ca 2+ dynamics (both Ca 2+ transients and Ca 2+ sparks) in shaping the force-frequency curve during the rapidly paced cardiomyocytes. This is an advance of previous work in common pool models, because the newly developed local control model simulated the Ca 2+ sparks governing excitationcontraction coupling dynamics [28]. The model displayed a positive force-frequency curve, starting from 0.2 Hz and ending at 4 Hz pacing; the calcium transients accompanied by increases to the RyR2 PO; RyR2 adaptation; and SR Ca 2+ load. The RyR2 Popen started to decline after 4 Hz pacing. Similarly, the upward trajectory of SR Ca 2+ ceased and started to decline slowly. However, the fraction of RyR2 channel was observed in the adapted state after 6 Hz. The positive FFR is the characteristic of the mammalian myocardium, including guinea pigs and rabbits [3,45]. Endoh [3] and Varian et al. [46] reported that, in rabbits, the Ca 2+ transient is positive until 4 Hz pacing and is negative at higher pacing frequencies. The model is used to dissect the molecular mechanism of FFR, as the processes governing the shape of the FFR are only moderately understood [47]. The model demonstrates that the interplay of the RyR2 Popen dependence on SR luminal Ca 2+ , SR Ca 2+ load, and RyR2 adaptation, which govern the shape of the FFR.
Initially RyR2 adaptation was thought to be relatively slow as a negative control mechanism, but Valdivia et al. [48] found it ~10-fold faster with Mg 2+ [49]. The adaptation rate for this model was 7 s −1 . Jafri et al. [28] reported that the adaptation of RyR2 decreased its Popen with the beginning of rapid pacing. In this model, the RyR2 adaptation increased as the pacing frequency increased, while RyR2 Popen showed a positive slope before 4 Hz and a negative slope above 4 Hz. In the model, we found 5% of RyR2 in adaptation state in 1 Hz; it reached 13% at 6 Hz during diastole, and the number increased

Discussion
The frequency-dependent performance of the myocardium changes the contractile strength of the heart. For humans, rabbits or guinea pigs, this force-frequency curve shows biphasic behavior, initially increasing and then decreasing [3,[5][6][7][8]. The force generated by the cardiac myocyte is a monotonically increasing function of the amplitude of the myoplasmic Ca 2+ transient [44]. A newly developed model of the ventricular myocytes of a guinea pig was used to study the role Ca 2+ dynamics (both Ca 2+ transients and Ca 2+ sparks) in shaping the force-frequency curve during the rapidly paced cardiomyocytes. This is an advance of previous work in common pool models, because the newly developed local control model simulated the Ca 2+ sparks governing excitation-contraction coupling dynamics [28]. The model displayed a positive force-frequency curve, starting from 0.2 Hz and ending at 4 Hz pacing; the calcium transients accompanied by increases to the RyR2 P O ; RyR2 adaptation; and SR Ca 2+ load. The RyR2 P open started to decline after 4 Hz pacing. Similarly, the upward trajectory of SR Ca 2+ ceased and started to decline slowly. However, the fraction of RyR2 channel was observed in the adapted state after 6 Hz. The positive FFR is the characteristic of the mammalian myocardium, including guinea pigs and rabbits [3,45]. Endoh [3] and Varian et al. [46] reported that, in rabbits, the Ca 2+ transient is positive until 4 Hz pacing and is negative at higher pacing frequencies. The model is used to dissect the molecular mechanism of FFR, as the processes governing the shape of the FFR are only moderately understood [47]. The model demonstrates that the interplay of the RyR2 P open dependence on SR luminal Ca 2+ , SR Ca 2+ load, and RyR2 adaptation, which govern the shape of the FFR.
Initially RyR2 adaptation was thought to be relatively slow as a negative control mechanism, but Valdivia et al. [48] found it~10-fold faster with Mg 2+ [49]. The adaptation rate for this model was 7 s −1 . Jafri et al. [28] reported that the adaptation of RyR2 decreased its P open with the beginning of rapid pacing. In this model, the RyR2 adaptation increased as the pacing frequency increased, while RyR2 P open showed a positive slope before 4 Hz and a negative slope above 4 Hz. In the model, we found 5% of RyR2 in adaptation state in 1 Hz; it reached 13% at 6 Hz during diastole, and the number increased during the systolic phase. This RyR2 in the adaptive state played the role in the lesser opening of RyR2 and the release of less SR Ca 2+ to the cytosol. After running simulations reducing adaptation by 20%, there was an increase in the Ca 2+ transient amplitude, the number of Ca 2+ sparks and their amplitudes. The RyR2 P open went up by 7% because both luminal dependency and stimulus by subspace Ca 2+ provide positive feedback to the P open . These results agree with the deterministic model developed by Jafri et al. [28], in that adaptation is important for producing FFR-related behavior. However, the peak of the FFR did not shift with the 20% reduction, suggesting that adaptation is not entirely responsible for the negative FFR. In the whole heart, other mechanisms have been suggested to be important. For example, Puglisi et al. [50] stated that, in the heart, it would not only be RyR2 adaptability, but the heart adapting to generate more force by increasing Ca 2+ load in the intracellular compartments at a high-frequency rate and increasing myofilament sensitivity [44]. For instance, in β-adrenergic stimulation, it prepares the ventricles to accommodate the higher beating rate.
It has been reported previously that the faster pacing rate leads to a physiological shortening of APD. Szigligeti et al. [45] changed positive FFR to negative FFR by shortening APD. The model displays this APD restitution with increased pacing frequencies ( Figure 4F). In other experiments, as the APD decreased with increased pacing frequency, the maximum rate of ventricular pressure was reciprocal to the APD (Franz, 1983). The frequency dependence of APD is caused by a decrease of the inward current, L-type current, and an increase of outward current, Na + -Ca 2+ exchange current [51]. The model also displayed a frequency-dependent decrease of L-type current (I LCC ) due to Ca 2+ -dependent inactivation and the increase of I ncx due to activation by Ca 2+ (Figure 5). The reduction in L-type current means that there are fewer triggering events for Ca 2+ sparks, which would reduce RyR2 P open as pacing frequency increases. This also contributes to the negative FFR. The local control model captures this mechanism that the previous common pool model could not.
The frequency-dependent increase of SR Ca 2+ load is the major contributor to positive slope FFR [3]. In guinea pigs, the SR Ca 2+ load increases with increased pacing rate. However, in rats, a negatively sloped FFR, the SR Ca 2+ load does not increase significantly with increased pacing frequency [52]. The model displayed increased SR Ca 2+ load with each pacing rate increase. The higher SR load produced larger Ca 2+ transients via upward trending RyR2 P open -modulated luminal dependency. This role of SR Ca 2+ has been illustrated previously by the positive staircase phenomenon [17]. An increase in heart rate increases the force of contraction generated by the myocyte; this phenomenon is associated with intracellular Ca 2+ handling in the myoplasm. In steady-state, with every depolarization, the influx of Ca 2+ from L-type channels leads to Ca 2+ release from SR. Myocyte relaxes when Ca 2+ returns to its original concentration by removing Ca 2+ from cytosol refill back to the SR by SERCA and efflux via I ncx . However, when the pacing frequency increases, the time averaged flux through L-type channels increases because there are more action potentials per second. The increased [Ca 2+ ] myo activates SERCA, which increases the sequestration of Ca 2+ into the SR. The increased SR increases the RyR2 P open in the model, similar to previous experiments [53]. Sobie et al. [54] reported RyR2 activity linearly depends on luminal Ca 2+ ; lowering the luminal value shifted the luminal regulation away from the RyR2, and the CICR-related activities were affected and a smaller number of Ca 2+ sparks, with smaller average amplitudes, were detected. The simulations shown in Figure 10 show that, when Ca 2+ efflux via the NCX is increased or Ca 2+ influx via the L-type Ca 2+ channel is decreased, thereby lowering the total myocyte Ca 2+ , there is limited filling of the SR, with an increase in pacing frequency resulting in a negative FFR.
A positive FFR is an intrinsic contractile property of a ventricular myocyte in mammals, and it is the result of a frequency-dependent acceleration of relaxation [3,5]. The model reproduced this behavior. Varian et al. [46] performed an FFR experiment in rabbits (in vivo) and found similar FFR behaviors as those reported in this article. They have shown that both FFR and Ca 2+ transient displayed a positive slope until they reached 4 Hz, and the current model followed a similar pattern. They did not report data at higher pacing frequencies, because it was thought that, with the high metabolic demand and greater rundown, data might be compromised. Endoh [3] also reported similar results from rabbit papillary muscle: the Ca 2+ transient was still positive until 4 Hz, but the contractile force associated with the amplitude of Ca 2+ transients showed positive FFR from 0.13 Hz to 3.30 Hz, then started to dissociate just before 4 Hz. He also observed negative FFR in rabbits at a higher frequency, and believed this occurred because of altered Ca 2+ handling and Ca 2+ overload. The model suggests that the SR loading with increased pacing rate leads to this positive slope of the FFR.
The increase in SR Ca 2+ load also increases contraction force generated by the myocyte in rapid pacing. The first Ca 2+ transient from lower to higher pacing becomes shorter, due to smaller recovery time to RyR2 from the previous inactivation, with a shortened diastolic phase. Continuous pacing at increasing frequencies leads to a gradual increase in a positive staircase before reaching a steady state. If the pacing rate decreases from a higher to lower frequency, opposite to the previous condition, the peak of first Ca 2+ transient becomes greater before gradually a falling to the steady state level. This is because of increased SR load due to an increased influx of Ca 2+ per unit time in rapid pacing, as well as enough time to reactivate RyR2 channels due to the elongated diastolic phase. We have also found three factors that contribute to a positive staircase, similar to Bers [17], in higher pacing: (a) increased L-type Ca 2+ current per unit time (not per unit beat); (b) higher diastolic [Ca 2+ ] myo ; and (c) increased SR Ca 2+ load [Ca 2+ ] sr available to be released in subsequent beats.
The FFR is an important indicator in finding failing or non-failing hearts. The twitch tension (FFR) rises in a normal heart, but does not rise in a failing one [55]. The heart failure is characterized by the decay of contraction and small systolic Ca 2+ transients [56]. Eisner et al. [56] also explained the small Ca 2+ transients in two theories: decreased activities of SERCA2a to reload SR and decreased RyR2 P open. Although SERCA2a activity was not the direct focus of the current research, with the negative FFR (after 5 Hz), we found decreased SR Ca 2+ load and reduced RyR2 P open , which both played a role in releasing the SR Ca 2+ in the smaller transients. The model also showed a reduction in I LCC and improved I ncx ; the former is responsible for decreasing RyR2 activation and the latter competes with SERCA2a to extrude more Ca 2+ out of the myoplasm. In the model, with the 20% dephosphorylation of the RyR2, we found a small decrease (5%) in the SR Ca 2+ spark leak, which increased SR Ca 2+ . In the failing heart, the opposite occurs, there is an increase in RyR2 P open due to hyperphosphorylation, but it also increases diastolic Ca 2+ leak and limits the amount of Ca 2+ in the SR [57]. In agreement with Endoh [3], the SERCA2a pumps in the failing heart exhaust their capacity to reload SR Ca 2+, and the positive FFR turns into negative FFR and the function of the heart is ceased [3]. The model also gave more reasons to believe that decreased SR load caused negative FFR, rather than SR leak. For a normal myocyte, the FFR plot indicates that an increase in pacing rate results in higher Ca 2+ levels in the myocytes and increases in the contractile force, but a dissociation in that force is necessary to prevent the myocyte from any mechanical damage. A healthy myocardium needs the positive FFR to continue a contractile behavior; in heart failure, the heart loses its ability to refill SR to continue to the frequency-dependent positive FFR [58].
RyR2 phosphorylation occurs during beta adrenergic stimulation, which can occur during normal physiology, such as exercise, or during diseases such as heart failure. The phosphorylation of RyR2 increases the diastolic Ca 2+ leak [57], and a decrease in RyR2 P open should decrease Ca 2+ leak. In comparing the model data of 1 Hz pacing, we found that the number of Ca 2+ sparks increase (1441 ± 914) with the reduction in RyR2 opening rate constant. The total diastolic leak rate for 1 Hz pacing was 28,687 ± 1632 for the control simulation, while it was 27,246 ± 1589 for the 20% reduced constant rate simulation. A 20% reduction in the RyR2 channel association rate (k + a ) for the transition from closed to open states showed that it has a negative effect on RyR2 P open ; however, in all other simulations, it always remains a fixed value, so it has no role in affecting the outcome. On the other hand, the luminal Ca 2+ and RyR2 adaptation play a positive and negative role in controlling RyR2 P open, respectively. However, in end-stage heart failure, the expression of SERCA is reduced and the expression of NCX is increased, with the consequence that the SR does not increase its Ca 2+ content during rapid pacing [59]. In the model, Figure 10 shows how increasing NSR Ca 2+ efflux abolishes the positive part of the FFR, matching the negative FFR [60] seen in heart failure [3].
The simulations in Figure 11 show that, when there is higher availability of SR Ca 2+ to be released, the open fraction of RyR2 is also high; conversely, if SR Ca 2+ load is lesser, the open fraction of RyR2 is also reduced. The interplay of these two performs a major role in the force-frequency relationship, and this is a widely studied phenomenon in both experimental and model settings. This protocol pacing is well studied in both experimental [17] and model settings [28], and the current model was also able to produce similar results. Due to the stochastic opening and closing of Ca 2+ channels, it was interesting to see the variations in the amplitude of RyR2 open fraction in the consecutive beats, even in a steady state. An increase in the open fraction rate was seen at the beginning of a lower frequency of pacing.

Conclusions
Simulations with the newly developed model for excitation-contraction coupling in guinea pig ventricular myocyte explored mechanisms behind the experimentally observed force-frequency relations. The continuous refilling of the SR by SERCA pump and enhanced SR Ca 2+ release is highly critical for increased force-frequency response. The model presented here predicts both the cellular and subcellular mechanism of FFR. In the cellular mechanism, the model suggests, in agreement with previous modeling studies by Jafri and co-workers [28], that diastolic sarcoplasmic reticulum (SR) [Ca 2+ ] SR and RyR2 adaptation increases with increased stimulation frequency, giving rise to increasing, rather than falling, amplitude of the myoplasmic [Ca 2+ ] myo transients. However, the model also suggests that the reduction of L-type Ca 2+ current with an increased pacing frequency due to Ca 2+dependent inactivation also contributes to the negative part of the FFR. The model was able to give such insight because it also allowed the dissection of these frequency-dependent mechanisms down to the spark level, gaining a deeper perspective into the regulation of cardiac calcium dynamics.

A.1. Calcium Release
The differential equations describing Ca 2+ at the release site are shown below, the details can be found in [27,28] where (i) is an index indicating the i tch specific Ca 2+ release unit, λ ds = V ds /V myo is the volume fraction that scale the fluxes, defined based on myoplasmic volume, to the subspace volume compartment. The membrane buffers used here are similar to those used previously [27][28][29].

A.2. Ryanodine Receptor Type-2 Model
This RyR2 model consists of 20,000 Ca 2+ release units in each cardiac cell and their state behavior is independent. Each Ca 2+ release unit contains a cluster of 49 stochastically gated RyR2s and 14 L-type channels. In the model, each Ca 2+ activated RyR2 is represented by two states-closed and open (see in Figure 1 Coupling of the RyR2 channels in the cluster uses allosteric coupling energies: ε OO , ε CC as before [27,28]. There is a possibility that after some time the channels might transition to an adapted state (C3) with the transition rare k + b . At each release site, the transition between two cluster states X i and X j is represented as where N O→C and N C→O represents the number of channels in cluster state X i and X j that can be switched from Open to Closed (Closed to Open). The coupled gating in the cluster is computed by the formula below with a j is the average allosteric connectivity for each RyR based on mean-field approach [61]. The variables N closed and N open represent the current number of non-conducting, and conducting channels in the cluster state X i , respectively. The non-junctional or 'rogue' RyR2 is modelled deterministically using the formula k + ryr−nj has the same functional form as k + ryr , except the subspace calcium [Ca 2+ ] ds in the formula of is now replaced by [Ca 2+ ] myo . In the model, the fraction of non-junctional RyR2s (ryr-nj) is 5%.

A.3. L-Type Ca 2+ Channel Model
The 6-state L-type Ca 2+ channel (LCC) model, Figure 2, is derived from 5-state LCC model for rat ventricle myocyte from Sun and co-workers with parameters adjusted for the new spark model [35]. Permeation of the L-type Ca 2+ channel used the Goldman-Hodgkin-Katz (GHK) formalism [60] with (i) represents the index of a release site N

A.4. Na + Channel Model
The whole-cell Na + current was derived from the formula with the ODE for gating variable follows the first-order differential equation where x Na represents the dimensionless gating variables, either m, h, or j, in the range between 0 and 1.0 and x ∞ is the steady-state value for the gating variable. The unit for the time constants is in seconds. There are five different K + channel currents (I ktof , I ktos , I K1 , I Kp and I Kss) . The formula for fast and slow transient outward currents (I Ktof , I Ktos ) respectively are based on the model developed using data observed in mouse [61]. with g K1 = 0.50; g Kss = 0.0421; g Kto f = 0.1598; g Ktos = 0.00669; g bK = 0.
with A m is the SL surface area and V myo is the myoplasmic volume The NCX is given by the formula with I ncx is the maximal NCX current. The parameters are given in Table A2.

A.7. Background Currents
We have background currents on the model were computed based upon the Ohm-law Three background currents are calcium, sodium, and potassium I bNa = g bNa (V m − E Na ) (A43) with the reversal potentials are derived from Nernst equation.
with z x is the valance of ion X.

A.9. Sarcoplasmic Reticulum Ion Pumps
The 2-state reduced model of the Tran-Crampin model [62] was chosen with parameters selected so that the maximum pumping rate per SERXA molecule, is 5 s −1 .

A.10. Buffering
The three endogenous buffers of calmodulin (CaM), troponin (Trpn), and the phospholipids of the SR membrane (SRm) are used for the bulk myoplasm. The buffering parameters are listed in Table A3.

A.11. Membrane Potential
Our AP model in guinea pig is based upon Luo and Rudy and Jafri and Rice [30,32]. In 1991 Luo Rudy published the first Luo-Rudy Model in ventricular myocytes of guinea pig [31].
During AP, the dynamics of the membrane voltage V m is governed by the equation I background = I bCa + I bNa + I bK (A58) NOTE: The units: membrane potential and reversal potential are in mV, while other parameters are given in the form of unit-density, e.g., the specific membrane capacitance C sc = 1 µF/cm 2 , the ionic currents are given in (µA/cm 2 ), except the L-type currents are measured at whole-cell level (µA) which needs to be converted to current density by dividing the membrane surface area A m (cm 2 ).