Effects of Sarcolemmal Background Ca2+ Entry and Sarcoplasmic Ca2+ Leak Currents on Electrophysiology and Ca2+ Transients in Human Ventricular Cardiomyocytes: A Computational Comparison

The intricate regulation of the compartmental Ca2+ concentrations in cardiomyocytes is critical for electrophysiology, excitation-contraction coupling, and other signaling pathways. Research into the complex signaling pathways is motivated by cardiac pathologies including arrhythmia and maladaptive myocyte remodeling, which result from Ca2+ dysregulation. Of interest to this investigation are two types of Ca2+ currents in cardiomyocytes: 1) background Ca2+ entry, i.e., Ca2+ transport across the sarcolemma from the extracellular space into the cytosol, and 2) Ca2+ leak from the sarcoplasmic reticulum (SR) across the SR membrane into the cytosol. Candidates for the ion channels underlying background Ca2+ entry and SR Ca2+ leak channels include members of the mechano-modulated transient receptor potential (TRP) family. We used a mathematical model of a human ventricular myocyte to analyze the individual contributions of background Ca2+ entry and SR Ca2+ leak to the modulation of Ca2+ transients and SR Ca2+ load at rest and during action potentials. Background Ca2+ entry exhibited a positive relationship with both [Ca2+]i and [Ca2+]SR. Modulating SR Ca2+ leak had opposite effects of background Ca2+ entry. Effects of SR Ca2+ leak on Ca2+ were particularly pronounced at lower pacing frequency. In contrast to the pronounced effects of background and leak Ca2+ currents on Ca2+ concentrations, the effects on cellular electrophysiology were marginal. Our studies provide quantitative insights into the differential modulation of compartmental Ca2+ concentrations by the background and leak Ca2+ currents. Furthermore, our studies support the hypothesis that TRP channels play a role in strain-modulation of cardiac contractility. In summary, our investigations shed light on the physiological effects of the background and leak Ca2+ currents and their contribution to the development of disease caused by Ca2+ dysregulation.


INTRODUCTION
Ca 2+ concentrations are dynamically controlled in cardiomyocytes by a complex regulatory system comprising ion channels, transporters, exchangers, regulatory proteins, and ion buffers. Intricately regulated levels of Ca 2+ concentrations are critical for electrical activity, excitationcontraction coupling, and other signaling pathways. In many cardiac diseases, the delicate balance of Ca 2+ cycling is perturbed. Ca 2+ dysregulation underlies maladaptive cardiac remodeling. A complete understanding of all components of Ca 2+ handling is essential for the development of therapeutic strategies to attenuate cardiac pathologies.
Many ion channels underlying the Ca 2+ signaling in cardiomyocytes are well characterized. Ca 2+ signaling related to excitation-contraction coupling primarily relies on sarcolemmal Ca 2+ entry through voltage-gated L-type Ca 2+ channels (LTCC) to trigger Ca 2+ -induced Ca 2+ release from the sarcoplasmic reticulum (SR) through ryanodine receptors (RyR), and, subsequently, Ca 2+ extrusion via the sodium-calcium exchanger (NCX) and reuptake into the SR through sarco/ endoplasmic Ca 2+ -ATPase (SERCA). Other Ca 2+ currents contribute to the modulation of Ca 2+ concentrations in cardiomyocytes as well. Of interest to this study are background Ca 2+ entry through the sarcolemma and SR Ca 2+ leak. Understanding of the physiological role of these Ca 2+ currents is still incomplete.
Beyond Ca 2+ transport across the sarcolemma from the extracellular space into the cytosol through LTCC and NCX, sarcolemmal Ca 2+ transport comprises background Ca 2+ entry. In resting cardiomyocytes, [Ca 2+ ] i is of the order of 100 nM, which indicates the existence of a background Ca 2+ entry pathway to balance the Ca 2+ efflux of NCX (Eisner et al., 2020). Resting ventricular myocytes depleted of Ca 2+ stores with caffeine are able to reload the SR by a mechanism that involves extracellular Ca 2+ , demonstrating further evidence of background Ca 2+ entry (Terracciano and Macleod, 1996). Based on studies measuring background Ca 2+ influx in rat ventricular myocytes of the order of 2-5 μmol/L per second (Choi et al., 2000) or 4 μmol/L per second (Sankaranarayanan et al., 2017), the background Ca 2+ entry is approximately 10% of the influx through LTCC current (5-10 μmol/L each action potential) at normal heart rates (Eisner et al., 2020). The identity of background Ca 2+ flux is still poorly defined, but several channels have been suggested as contributors.
One study identified a Ca 2+ entry mechanism that is blocked by the nonspecific agent gadolinium (Gd 3+ ) (Kupittayanant et al., 2006). Connexin hemichannels are candidates for background Ca 2+ entry since they can be inhibited by Gd 3+ (Stout et al., 2002). While connexin hemichannels primarily form pairs to allow ion fluxes between cells at intercalated discs, some are present as hemichannels in the surface membrane of a single cell (Wang et al., 2012;Leybaert et al., 2017) and may, therefore, provide a route for Ca 2+ entry. However, primary candidates for background Ca 2+ entry include members of the family of Transient Receptor Potential (TRP) channels, which are also sensitive to Gd 3+ . The mdx mouse model of muscular dystrophy exhibits [Ca 2+ ] i and [Na + ] i overload that can be blocked by Gd 3+ , and the increase in cation entry has been suggested to involve TRPC channels (Mijares et al., 2014). Myocytes from old mdx mice exhibit increased expression of a putative stretch-activated channel (SAC), TRPC1. Elevated [Ca 2+ ] i levels can also be reduced to [Ca 2+ ] i levels of WT myocytes when exposed to SAC blockers streptomycin or GsMTx-4 (Williams and Allen, 2007;Ward et al., 2008). Upregulated TRPC1 also contributes to increased [Ca 2+ ] i through SAC in hypertrophic myocardium of rats following isoproterenol injection (Chen et al., 2013). Background Ca 2+ entry involved in maladaptive cardiac remodeling was more recently shown to critically depend on both TRPC1 and TRPC4 (Camacho Londono et al., 2015, Camacho Londoño et al., 2021. TRPC6 channels have also been shown to modulate cytosolic Ca 2+ transients and SR Ca 2+ load through sarcolemmal Ca 2+ entry . Furthermore, TRPV4 can modulate Ca 2+ transients and SR load, and participates in hypoosmotic stress-induced cardiomyocyte Ca 2+ entry (Rubinstein et al., 2014;Jones et al., 2019).
Ca 2+ transport across the SR membrane from the intracellular store into the cytosol is known as SR Ca 2+ leak. During an action potential, activation of RyR clusters triggered by Ca 2+ entry through LTCC results in a synchronized release of a large amount of Ca 2+ from the SR that forms the Ca 2+ transient and leads to cardiomyocyte contraction. Activation of RyR clusters causes Ca 2+ sparks. These sparks are an important pathway for SR Ca 2+ leak. RyR Ca 2+ leak can occur also through a mechanism independent of sparks (Santiago et al., 2010). Further, total SR Ca 2+ leak includes a component separate from RyRs (Zima et al., 2010). Many candidates have been suggested as components of this leak, yet it is still ill-defined. A candidate is the inositol 1,4,5-trisphosphate (IP 3 ) receptor (IP 3 R), which is a Ca 2+ release channel expressed at lower densities than RyRs in cardiomyocytes and found upregulated in heart failure (HF) (Go et al., 1995;Ai et al., 2005) and therefore may be a relevant contributor to SR Ca 2+ leak (Zima et al., 2010). Members of the TRP family have also been suggested to contribute to SR Ca 2+ leak. TRPC1 was found to operate as a SR Ca 2+ leak channel in skeletal muscle (Berbey et al., 2009) and more recently in cardiomyocytes (Hu et al., 2020). Additional evidence suggests contribution of TRPC6, TRPM8, TRPP2, and TRPV1 to endoplasmic reticulum Ca 2+ leak in various cell types, although characterization is still incomplete for cardiomyocytes (Lemos et al., 2021).
TRP channels constitute primary candidates for explaining background and leak Ca 2+ currents through both the sarcolemma and the SR membrane. Interestingly, many members of the TRP family are known to be modulated by stretch (Inoue et al., 2009;Reed et al., 2014;Peyronnet et al., 2016). Cardiomyocyte contractility is known to respond to mechanical stretch in two phases: a rapid and a slow response (Calaghan and White, 1999). The rapid response is the cellular basis for the Frank-Starling Mechanism (FSM). It relies primarily on myofilament overlap and alteration of myofilament Ca 2+ sensitivity, and does not involve changes in Ca 2+ transients. The slow response has been termed slow force response (SFR) or stress-induced slow increase in contractility (SSC), describing the gradual increase in twitch force corresponding to an increase in [Ca 2+ ] i transients that develops over several minutes when stretch is sustained. Members of the TRP family are likely candidates for mechanotransduction of the SFR.
Though background and leak Ca 2+ currents are still poorly defined, they may play important roles in regulating Ca 2+ homeostasis and contractility and altered Ca 2+ dynamics in cardiac disease. The lack of complete understanding of the identity and mechanics of these channels, in addition to their relatively small amplitude compared to the voltage-gated Ca 2+ channels, make them difficult to study in vivo. A computational model provides the advantage to study the currents and their roles in isolation from other cellular mechanisms. In this study, we use a mathematical model of a human ventricular myocyte to analyze the individual contributions of background Ca 2+ entry and SR Ca 2+ leak to the modulation of Ca 2+ transients and SR Ca 2+ load. We also assess the effects of the Ca 2+ currents on cellular electrophysiology.

Mathematical Model of Ventricular Myocyte
We applied a mathematical model of a human ventricular myocyte (Grandi et al., 2010). The model and subsequent analyses were executed in MATLAB (R2020b). The model includes subsarcolemmal and junctional compartments beyond the cytosol compartment. Total background Ca 2+ current (I Cabk ) through the sarcolemma is defined as a summation of junctional (I Cabk, junc ) and subsarcolemmal (I Cabk,sl ) components: where F junc = 0.11 and F sl = 0.89 are constants that determine the fraction of total background current corresponding to the junctional and subsarcolemmal spaces, respectively, G bkg is the maximum conductance (5.513e-4 A/F) of the channels, V m is the transmembrane voltage, and E Ca,junc and E Ca,sl are Nernst potentials corresponding to the junctional and subsarcolemmal spaces, respectively. SR Ca 2+ leak (J leak ) describes the Ca 2+ flux out of the SR: where K leak is the leak constant (5.348e-6/ms), and [Ca 2+ ] SR and [Ca 2+ ] j describe the Ca 2+ concentrations in the SR and junctional space, respectively. In this study, we modulated G bkg and K leak independently and in conjunction to investigate the effects of altered sarcolemmal Ca 2+ entry and altered SR Ca 2+ leak, respectively. A graphical summary of the Ca 2+ currents in the cell model is shown in Figure 1. The total Ca 2+ current into the junctional (I Ca,tot,junc ) and subsarcolemmal (I Ca,tot,sl ) compartments are determined by Eqs 5, 6, respectively: I Ca tot sl I Ca sl + I Cabk sl + I pCa sl − 2 · I NCX sl , where I Ca describes the L-type Ca 2+ current, I pCa describes the sarcolemmal Ca 2+ pump current, and I NCX describes the NCX current.
FIGURE 1 | Diagram of Ca 2+ handling components in human ventricular myocyte model (Grandi et al., 2010). In this mathematical model, the cytosol is divided into junction, sub-sarcolemmal and bulk compartments. Ca 2+ currents through the sarcolemma enter both the junctional and sub-sarcolemmal compartments before diffusing to the bulk cytosol.
where C mem is the membrane capacitance, F is Faraday's constant, V junc is the volume of the junctional compartment, V sl is the volume of the subsarcolemmal compartment, V myo is the volume of the bulk cytosol, and V sr is the volume of the SR. J Ca,juncsl is the rate of Ca 2+ flux from junctional to subsarcolemmal compartments and J Ca,slmyo is the rate of Ca 2+ flux from the subsarcolemmal compartment to bulk cytosol. J RYR describes the SR Ca 2+ release and J SERCA describes the SR Ca 2+ pump. Ca 2+ buffering is described by J Ca,B,junction in the junctional compartment, J Ca,B,sl in the subsarcolemmal compartment, J Ca,B,cytosol in the bulk cytosol, and dCsqn b /dt in the SR. We refer to (Grandi et al., 2010) for more details and the equations describing the mathematical model.

Evaluating the Effects of Background and Leak Ca 2+ Currents on Ca 2+ Concentrations
To investigate the effects of altered background Ca 2+ entry, we modulated G bkg from 0 to 300% of its default value, 5.513e-4 A/F. To investigate the effects of altered SR Ca 2+ leak, we modulated K leak from 0 to 300% of its default value, 5.348e-6 ms −1 . Simulation durations were 1 min to establish steady state. We analyzed the last beat. We performed sensitivity analyses for the modulation of G bkg and K leak independently on measurements of resting V m , maximum V m , action potential duration to 90% repolarization ( where x is the fraction (%) of default G bkg or K leak . We evaluated the factor of the quadratic term, A, the linear term, B, and the constant term, C, of this fit for relative comparisons of the effects of modulated background Ca 2+ current or modulated SR Ca 2+ leak. This analysis was repeated for the model ran at 1, 2, and 3 Hz electrical excitation. We subsequently performed a dual-sensitivity analysis for the modulations of G bkg and K leak to understand how the two Ca 2+ currents interact and influence Ca 2+ handling in the cardiomyocyte together.

Qualitative Changes in V m , [Ca 2+ ] i , and [Ca 2+ ] SR When Background Ca 2+ Entry is Modulated
We first examined the changes in V m , [Ca 2+ ] i , and [Ca 2+ ] SR following modulation of the background Ca 2+ entry current, I Cabk, at pacing rates of 1, 2 and 3 Hz ( Figure 2). Changes to features of the action potential were minimal ( Figure 2A). Resting V m exhibited a positive relationship with % G bkg ( Figure 2B), while maximum depolarization exhibited a negative relationship ( Figure 2C) with increasing background Ca 2+ current. APD 90 showed a positive relationship, where increased background Ca 2+ entry lengthens the time for repolarization ( Figure 2D). We measured diastolic and systolic values, as well as amplitude of [Ca 2+ ] i transients ( Figure 2E). All measurements demonstrate positive relationships with increasing background Ca 2+ entry; the increase in systolic [Ca 2+ ] i ( Figure 2F) is greater than the increase in diastolic [Ca 2+ ] i ( Figure 2G), especially at the slowest pacing rate, so the amplitude increases a considerable amount ( Figure 2H). We measured the same parameters of [Ca 2+ ] SR transients, which also demonstrate an increase in amplitude caused by increased background Ca 2+ entry ( Figures 2I,L). The release was initiated from increased diastolic [Ca 2+ ] SR ( Figure 2J) and ended in reduced systolic [Ca 2+ ] SR ( Figure 2K). The relationships between background Ca 2+ entry and features of the [Ca 2+ ] SR transient were also augmented at the slowest pacing rate.

Qualitative Changes in V m , [Ca 2+ ] i , and [Ca 2+ ] SR When SR Ca 2+ Leak is Modulated
Changes following independent modulation of the SR Ca 2+ leak current, J leak, were also examined ( Figure 3). The effects of modulating K leak on V m were small ( Figure 3A). Resting V m exhibited a slight positive relationship ( Figure 3B). Maximum V m increased with increasing K leak at 1 Hz pacing but decreased with K leak at 2 and 3 Hz pacing ( Figure 3C). Repolarization time measured as APD 90 was unaltered by modulations of K leak ( Figure 3D). Diastolic [Ca 2+ ] i was also relatively unaffected by modulations of K leak (Figures 3E,F) (Figure 3K), both contributing to reduced [Ca 2+ ] SR amplitude with increased K leak ( Figure 3L).

Summary and Comparison of Independent Modulation of Background and Leak Currents
For comparison of the effects of altered G bkg or K leak on features of the action potential, we normalized the measurements of resting V m , maximum V m , and APD 90 to the default measurements from the model for the given pacing frequency ( Figures 4A-F). The relationships of measured vs. modified parameter, both represented as fraction (%) vs. default values, were fit to a 2nd order polynomial model for quantification of the effects (Figures 4G-I). The quadratic term of the fit is negligible for resting V m and maximum V m , demonstrating that these measurements exhibit primarily linear relationships with G bkg and K leak ( Figure 4G). However, the polynomial fit of APD 90 to G bkg has a strong quadratic term, demonstrating a non-linear relationship between APD 90 and G bkg ( Figure 4G). Since the relationships of resting V m and maximum V m are primarily linear, the linear term of the quadratic polynomial fit demonstrates the sensitivity of the measurements to changes in G bkg or K leak (Figure 4H). Both G bkg and K leak modulation have a positive relationship with resting V m , but the changes are negligible for K leak and minimal for G bkg , never exceeding an increase greater than 2% of the default model's resting V m . Maximum V m had a strong negative linear relationship with G bkg for all pacing frequencies, especially at 1 Hz pacing. Maximum V m had a positive linear relationship with K leak at 1 Hz pacing but negatives linear relationship for 2 and 3 Hz. The constant of the quadratic polynomial model represents the measurements for the modulated currents set to 0 ( Figure 4I). These results revealed that resting V m was largely unchanged by pacing rate. Maximum V m was slightly increased with no I Cabk and slightly reduced with no J leak at 1 Hz and slightly increased with no J leak at 3 Hz. APD 90 decreased without I Cabk but remained unchanged without K leak . Figure 5 contains a summary of [Ca 2+ ] i normalized by measurements of the default model for the sensitivity analyses of G bkg and K leak modulation. The relationships of measured parameter vs. modified parameter, both represented as % default values, were fit to a 2nd order polynomial model for quantification of the effects. The quadratic term indicates nonlinearity of the relationship ( Figure 5G). The quadratic term of the diastolic [Ca 2+ ] i fits were marginal, indicating a primarily linear relationship. For both G bkg and K leak modulations, we noticed a nonlinearity associated with systolic [Ca 2+ ] i and consequently the [Ca 2+ ] i amplitude, with the greatest degree of nonlinearity associated with the lowest pacing frequency. The linear term of the quadratic fit showed a strong linear sensitivity of the measured parameter to changes in G bkg or K leak ( Figure 5H). The positive sign of this term for G bkg modulations for each measurement, diastolic, systolic, and amplitude, demonstrates the positive relationship of these parameters with G bkg , with greater slope of the relationship for systolic and amplitude measurements.  Figure 6. Again, a 2nd order polynomial model for the normalized sensitivity analyses provided information about the relationships. The quadratic term of the fits demonstrates that systolic [Ca 2+ ] SR and [Ca 2+ ] SR transient amplitude both exhibited nonlinearity in the relationship to modulated G bkg ( Figure 6G). The linear term of the relationships demonstrated the best representation sensitivity of the measurements to G bkg or K leak modulation ( Figure 6H). The sign of this term was opposite for each measured parameter for G bkg vs. K leak modulation but similar in amplitude. With a thorough understanding of how G bkg and K leak modulations independently affect features of the action potential, [Ca 2+ ] i and [Ca 2+ ] SR transients, we subsequently performed a dual-parameter sensitivity analysis of G bkg and K leak together for 1 Hz pacing (Figure 7). Resting V m is affected minimally by G bkg , and negligibly by K leak , apparent by the nearly horizontal contour lines ( Figure 7A). Both increasing G bkg and increasing K leak resulted in increased V m , but increasing G bkg makes a larger contribution. Maximum V m was modulated in opposite directions by G bkg and K leak (Figure 7B). Increasing G bkg reduces maximum V m while increasing K leak increased maximum V m , but G bkg is the slightly more dominant effect. APD 90 experienced the greatest increase at large increases in G bkg and small K leak ( Figure 7C). Diastolic [Ca 2+ ] i was primarily affected by a positive relationship to G bkg , while K leak appeared to make no contribution ( Figure 7D) (Figure 8).

DISCUSSION
In this study, we evaluated the contributions of background Ca 2+ entry and SR Ca 2+ leak to V m and Ca 2+ concentrations in the SR and cytosol in cardiomyocytes in silico. Our investigations shed light on the differential effects of background and leak Ca 2+ currents in physiology, and also provide insight into their contributions to disease development due to Ca 2+ dysfunction. Below, we discussed background Ca 2+ entry as a mechanism to positively modulate Ca 2+ entry and SR Ca 2+ leak as a critical balancing mechanism to maintain homeostasis.

Background Ca 2+ Entry Positively Modulates Ca 2+ Concentrations
Our results show that background Ca 2+ entry has a positive relationship with diastolic [Ca 2+ ] i and [Ca 2+ ] SR , and the amplitude of their transients (Figure 2). The small increase in diastolic [Ca 2+ ] i is the direct result of an increase in background Ca 2+ entry, and the increase in [Ca 2+ ] SR follows as SERCA responds to pump the extra Ca 2+ into the SR. The small increase in diastolic [Ca 2+ ] i amplifies Ca 2+ -induced-Ca 2+ release through RyRs, explaining the increased transient amplitudes. This supports the physiological role of background Ca 2+ entry in increasing Ca 2+ concentrations. These results replicate prior findings for TRP family members suggested as Ca 2+ entry channels. For example, we provided evidence for TRPC6 as background Ca 2+ entry in neonatal rat ventricular   (Seo et al., 2014;Yamaguchi et al., 2017;Yamaguchi et al., 2018). These studies highlight one potential role for Ca 2+ entry in strained myocytes. Many of the candidates suggested as background Ca 2+ entry channels are known to be modulated by stretch (Inoue et al., 2009;Reed et al., 2014;Peyronnet et al., 2016), so strained myocytes would exhibit an increase in background Ca 2+ entry, leading to elevated diastolic levels and larger transients. This mechanism likely contributes to the SFR, increased contractile force following sustained stretch.  Figure 8 demonstrate that modulations in G bkg and K leak can be balanced in a way to cancel out the large opposing effects they have on Ca 2+ transient amplitudes ( Figure 8). This balancing mechanism of SR Ca 2+ leak could be critical to prevent Ca 2+ overload in the cell. Leak in the form of RyR sparks has been demonstrated as SR load regulator to prevent overload, with a steep dependency on [Ca 2+ ] SR (Shannon et al., 2002). However, large Ca 2+ release events through RyR sparks increase sensitivity for arrhythmia (George, 2008). RyR sparks are most prevalent at high [Ca 2+ ] SR , but nonspark RyR leak and non-RyR leak do not appear to exhibit the same steep dependency on [Ca 2+ ] SR and therefore might function differently. Thus, a mechanism of non-RyR leak may be to regulate compartmental Ca 2+ before the cells become overloaded, and RyR sparks increase as a more extreme measure. Like channels involved in background Ca 2+ entry in cardiomyocytes, candidates for SR Ca 2+ leak also include members of the TRP family and were suggested to be mechanomodulated (Inoue et al., 2009;Reed et al., 2014;Peyronnet et al., 2016). This indicates that the background and leak Ca 2+ currents could be modulated in conjunction. Non-RyR leak that can be modulated by stretch may provide a more moderate and steady regulation on a beat-by-beat basis in conjunction with background Ca 2+ entry modulated by stretch. Recently we demonstrated that TRPC1 constitutes an SR Ca 2+ leak channel, and its overexpression resulted in decreased SR Ca 2+ load (Hu et al., 2020). TRPC1 channels are suggested to be modulated by stretch, indicating that the reduction in SR Ca 2+ load could be a regulatory mechanism to match increased background Ca 2+ entry through, e.g., TRPC6 channels. We speculate that background Ca 2+ entry and SR Ca 2+ leak fulfill a critical homeostatic function in the modulation of Ca 2+ concentrations throughout the cardiomyocyte in response to strain.

Background and Leak Ca 2+ Currents May Contribute to Hypertrophy and HF Under Chronic Pressure Overload
Both background Ca 2+ entry and SR Ca 2+ leak through TRP channels are likely to be modulated by cardiomyocyte strain (Inoue et al., 2009;Reed et al., 2014;Peyronnet et al., 2016). Under chronic pressure overload conditions, strain-modulation of TRPC channels could increase background Ca 2+ entry and SR Ca 2+ leak, and thus dysregulate Ca 2+ . Cardiac disease is perpetuated by Ca 2+ dysregulation, and a stray from its homeostatic balance. Some of the suggested ion channels for these Ca 2+ currents were found to be upregulated in models of cardiac disease, suggesting a role in pathogenesis (Ahmad et al., 2017;Hof et al., 2019 (Bers, 2006). Based on the demonstration of a balancing mechanism between background Ca 2+ entry and SR Ca 2+ leak in this study, it is reasonable to speculate a difference between maintenance of this balance in HFpEF vs. a stray from this balanced relationship towards overcompensated leak in HFrEF. In addition to reducing SR Ca 2+ available for release, causing systolic dysfunction, increased SR Ca 2+ leak can be problematic, e.g., triggering arrhythmias and being energetically costly due to increased use of ATP to repump Ca 2+ (Bers, 2014 (Frampton et al., 1991;Antoons et al., 2002;Dibb et al., 2007;Horváth et al., 2017;Sankaranarayanan et al., 2017). Background Ca 2+ entry and SR Ca 2+ leak also both exhibit a frequency effect. The parameters we measured are all more sensitive to modulations of the Ca 2+ currents at slower pacing rates than at faster pacing rates (Figures 5, 6). The sensitivity of each measured [Ca 2+ ] i and [Ca 2+ ] SR parameter to G bkg and K leak is greatest in amplitude for 1 Hz pacing. An explanation is that at slower pacing rates, the background and leak currents have relatively more time to contribute to the total Ca 2+ flux per beat vs. the voltage-gated ion channels that open during the action potential and are closed at rest.

Modulation of I Cabk and J leak has Marginal Effects on Action Potentials
Modulating K leak had negligible effects on the action potential for any pacing frequency (Figure 4). For the values of K leak tested, we found that the SR Ca 2+ leak flux does not significantly contribute to sarcolemmal electrophysiology. Modulating G bkg has marginal effects on features of action potentials (Figure 4). While G bkg positively correlates with increased resting V m , an increase to 300% G bkg only resulted in <2% change from basal resting V m . This minimal change in resting potential is unlikely to be functionally relevant. An increase to 300% G bkg also reduces maximum depolarization by 7% for 1 Hz pacing. The largest effect is an increase in action potential duration (APD90) by around 10% for maximal G bkg modulation. It has also been shown that APD prolongation leads to increased Ca 2+ (Bouchard et al., 1995), suggesting a positive feedback loop for electrical and Ca 2+ signaling. Another important note is that APD increase is known to be inotropic, e.g., in rat ventricular myocytes (Bouchard et al., 1995). This indicates another mechanism for the contribution of background Ca 2+ entry to contractility. Conversely, prolonged APD can induce torsades de pointes tachycardia, leading to lifethreatening ventricular fibrillation (Roden and Hoffman, 1985;Ravens and Cerbai, 2008).

Limitations
Mathematical modeling of cellular electrophysiology provides a valuable resource for studying how aspects of cellular physiology interact and affect one another. It provides a means to investigate questions that cannot be easily answered in vivo. However, there are also caveats of mathematical modeling that should be considered. It should be noted that the definitions of I Cabk and J leak used in this model are general simplifications and meant to reproduce poorly defined currents. The equations lack specific gating conditions of the currents. The current equations were not parameterized to match experimental data which is only incompletely characterized in human ventricular myocytes. Instead, the current equations are adjusted such that the model reproduces overall physiological action potentials and calcium transients. This is an important consideration, since the magnitudes of these currents could be largely different in living cells. Thus, interpreting the results of this study should focus on the qualitative trends. As the specific ion channels that contribute to Ca 2+ entry and leak are identified and characterized, future work can aim to refine the current definitions and provide detailed current models to replace the general simplifications of I Cabk and J leak .
In a similar way, other ion currents in the model are not fully defined. For example, some K + channels and isoforms of the Na + / K + -ATPase are modulated by localized Ca 2+ concentrations (Tian and Xie, 2008;Weisbrod, 2020). However, the model does not include any Ca 2+ -dependent terms in the definitions of these currents. The inclusion of these interactions may alter the effects we see on V m in this study. Future work could address this limitation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
MS and FS designed the study. MS implemented the modeling, analyzed simulation data and drafted the manuscript. All authors critically revised the manuscript and approved the version to be published.

FUNDING
This work was supported by Nora Eccles Treadwell Foundation.