Cycling of Ca2+ between the sarcoplasmic reticulum (SR) and myoplasm is an important component of skeletal muscle resting metabolism. As part of this cycle, Ca2+ leaks from the SR into the myoplasm and is pumped back into the SR using ATP, which leads to the consumption of O2 and generation of heat. Ca2+ may leak through release channels or ryanodine receptors (RYRs). RYR Ca2+ leak can be monitored in a skinned fiber preparation in which leaked Ca2+ is pumped into the t-system and measured with a fluorescent dye. However, accurate quantification faces a number of hurdles. To overcome them, we developed a mathematical model of Ca2+ movement in these preparations. The model incorporated Ca2+ pumps that move Ca2+ from the myoplasm to the SR and from the junctional space (JS) to the t-system, Ca2+ buffering by EGTA in the JS and myoplasm and by buffers in the SR, and Ca2+ leaks from the SR into the JS and myoplasm and from the t-system into the myoplasm. The model accurately simulated Ca2+ uptake into the t-system, the relationship between myoplasmic [Ca2+] and steady-state t-system [Ca2+], and the effect of blocking RYR Ca2+ leak on t-system Ca2+ uptake. The magnitude of the leak through the RYRs would contribute ∼5% of the resting heat production of human muscle. In normal resting fibers, RYR Ca2+ leak makes a small contribution to resting metabolism. RYR-focused pathologies have the potential to increase RYR Ca2+ leak and the RYR leak component of resting metabolism.
Cycling of calcium ions (Ca2+) between the sarcoplasmic reticulum (SR) and the myoplasm is central to excitation–contraction coupling in muscle: release of Ca2+ into the myoplasm and its subsequent binding to sites on troponin C (TnC) initiate and maintain contraction, and removal of Ca2+ from the myoplasm initiates relaxation. However, Ca2+ cycling is important even in unstimulated or resting muscle and contributes to resting metabolism. In resting muscle, there appears to be a steady leak of Ca2+ from the SR into the myoplasm that is balanced by removal of Ca2+ from the myoplasm into the SR via the SR Ca2+ pump (Chinet et al., 1992). The concentration of free, or unbound, Ca2+ in the myoplasm is steady in resting muscle (Konishi and Watanabe, 1995), indicating that the rate at which Ca2+ leaks into the myoplasm is equal to the rate at which it is pumped into the SR. The SR Ca2+ pump (SERCA) is powered by ATP and, in a metabolic steady state, ATP is replenished by mitochondrial oxidative phosphorylation at the same rate at which it is used. The ATP regenerating processes produce heat, and the heat associated with SERCA activity accounts for about one-quarter of the heat produced by resting mouse skeletal muscles (Chinet et al., 1992).
Two Ca2+ leak processes, that are not mutually exclusive, have been proposed. First, Ca2+ can leak from the SR through the Ca2+ release channels or ryanodine receptors (RYRs; Cully et al., 2016). The importance of this pathway is that the leakiness of RYRs has been shown to vary across a range of naturally occurring variants of RYRs, perhaps most strikingly in the highly modified heater muscles of some fish (Morrissette et al., 2003), and to be associated with pathological conditions, such as malignant hyperthermia (Cully et al., 2018). Second, Ca2+ can leak from the SR into the myoplasm via a mechanism independent of RYRs but that involves SERCA (Simonides and Van Hardeveld, 1988; Macdonald and Stephenson, 2001; Lamboley et al., 2014). The focus of the current study is the first of these mechanisms, Ca2+ leak through the RYRs.
Launikonis and colleagues have developed a method to measure the rate of Ca2+ leak through the RYRs in unstimulated or resting muscle fibers (Cully et al., 2016; Cully et al., 2017; Cully et al., 2018). A skinned muscle fiber preparation is used in which at least some of the Ca2+ that leaks through the RYRs is pumped into the t-tubules, which are sealed from the extracellular space in this particular preparation and are loaded with a Ca2+-sensitive dye. The inside of the fiber preparation is heavily buffered against changes in [Ca2+] by the presence of EGTA. In the junctional space (JS) adjacent to the RYRs, the presence of EGTA attenuates changes in Ca2+ concentration due to RYR leaks so that the variation in the concentration of unbound Ca2+ in the JS (CaJS) due to leaks is within the working range of the Ca2+ pump in the neighboring t-tubule membrane. In this way, the rate at which the pump (a plasma membrane Ca2+ ATPase [PMCA]) moves Ca2+ into the t-tubule reflects CaJS and thus the rate at which Ca2+ leaks through the RYRs. If Ca2+ leaks through the RYRs rapidly enough to raise CaJS, then the rate and extent of Ca2+ accumulation in the t-system will be increased. A key aspect of these experiments is that the RYR Ca2+ flux can be inhibited pharmacologically, using tetracaine. The magnitude of trans-RYR Ca2+ leak can then be inferred from the difference in the rate and extent of t-tubule Ca2+ accumulation before and after the inhibition (Cully et al., 2018). Those experiments indicate that there is a continual leak of Ca2+ through the RYR in unstimulated fibers of human skeletal muscle. However, the complexity of the experimental preparation makes it difficult to determine the magnitude of the RYR Ca2+ leak.
Specific aspects of the preparation and the experimental methods that complicate the relationship between RYR Ca2+ leak and t-tubule Ca2+ uptake are as follows: (1) the skinned fiber contains a high concentration of EGTA, which is used to control the myoplasmic Ca2+ concentration (CaM) to match the range over which the PMCA functions. However, the buffering action of EGTA must also attenuate the increase in CaJS that occurs when Ca2+ leaks into the JS. (2) It is likely that there is a leak of Ca2+ from the t-system into the myoplasm (Cully et al., 2018; Meizoso-Huesca and Launikonis, 2021). In that case, the empirical relationship between the concentration of unbound Ca2+ in the t-system (CaT) and CaJS reflects the combined characteristics of Ca2+ uptake by the PMCA and Ca2+ leak from the t-system. (3) JS volume (VJS) is small compared to that of the SR (VSR; VJS < 0.1% VSR; Eisenberg, 1983), so a leak of Ca2+ through the RYRs that makes only a small change in CaSR will produce a large change in the total Ca2+ concentration (i.e., free + bound to EGTA) in the JS. (4) The experimental protocol starts with emptying Ca2+ from the t-tubules and SR. It is not known whether the subsequent time course of SR filling, which in turn is likely to affect RYR leak, influences the time course of t-system Ca2+ accumulation.
One way of assessing how these aspects of the preparation and protocol influence the experimental results is to use a mathematical model of the preparation. The purpose of the current study was to develop a model that could be used to quantify the leak of Ca2+ from the SR through the RYRs, to identify factors that influence t-tubule Ca2+ accumulation, and to estimate the contribution of trans-RYR Ca2+ flux to the metabolism of resting muscle.
Materials and methods
Overview of the model
The model describes the movements of Ca2+ among the SR, t-tubule, JS, and myoplasm and is illustrated in Fig. 1. The subcellular system that was modeled included seven compartments that Ca2+ can occupy, four of which are physically discrete spaces and three of which are molecular compartments (i.e., Ca2+ in equilibrium with buffer molecules). The compartments are (1) t-tubules, (2) SR, (3) JS (i.e., region between t-tubule and SR Ca2+ release channels), (4) myoplasm, (5) EGTA in the myoplasm, (6) EGTA in the JS, and (7) Ca2+ buffers, including calsequestrin (Csq), in the SR.
Mathematical description of the model
The rates of change in concentration of Ca2+ in each of the compartments, taking into account the constraints on Ca2+ movements described in the legend of Fig. 1 and the relative volumes of the physical compartments, were described by ordinary differential equations. The equations incorporated Ca2+ movements due to (1) ATP-dependent pumping of Ca2+ from the myoplasm into the SR via SERCA, (2) ATP-dependent pumping of Ca2+ from the JS into the t-tubules via PMCA, (3) exchange of Na+ and Ca2+ between the JS and t-tubule via the sodium/calcium exchanger (NCX), and (4) diffusive exchange of Ca2+ between the JS and myoplasm. The concentrations of unbound Ca2+ in the myoplasm and JS were determined from the equilibrium between Ca2+ and EGTA, and that in the SR from the equilibrium between Ca2+ and the SR buffers. The equation formulation ensured that total Ca2+ in the system remained constant.
To provide the model with the complexity and responsiveness to local conditions that is evident in physiological systems, where possible mechanistic descriptions of model components were used. Models for SERCA and PMCA were taken from models described and validated in the literature. Cully et al. (2018) showed that t-system Ca2+ accumulation in human fibers was unaffected by pharmacological inhibition of NCX. Therefore, NCX was not incorporated into the final model. Temperature-sensitive model parameters obtained from the literature were adjusted to values appropriate for a temperature of 25°C, similar to that used for the muscle fiber experiments, assuming a Q10 (change in rate for a 10°C change in temperature) of 2.
SR Ca2+ pump (SERCA)
In each pump cycle, SERCA moves two Ca2+ from the myoplasm up a concentration gradient into the SR at the expense of one ATP (Weber et al., 1966; Makinose and Hasselbach, 1971). The pump was modeled as described by Tran et al. (2009). Those authors derived a 2-state model (Fig. 2 A) from a 12-state biochemical model. In reducing the model to 2 states, the relevant complexity of the 12-state model was maintained in that the apparent rate constants describing the forward and reverse transitions between the 2 states depend on pH, the concentrations of ATP and its hydrolysis products, and [Ca2+] on both sides of the SR membrane (see Appendix 1). The model is thermodynamically constrained; that is, the rate of Ca2+ pumping depends on the difference in free energy produced by ATP hydrolysis (ΔGATP) and the free energy generated by the Ca2+ concentration gradient (ΔGCa). The equations describing the pumping rate are given in Appendix 1.
To align the SERCA model with human skeletal muscle, appropriate values for the pump’s concentration and Ca2+ sensitivity were determined. Sensitivity of the pump to Ca2+ is quantified by its Ca50 value (the CaM at which the pump rate is 50% of its maximum). A constraint on the sensitivity is that the rate of Ca2+ pumping when CaM is at resting levels should not result in SERCA producing heat from oxidative regeneration of ATP consumed by the pump at a rate greater than the experimentally determined Ca2+-related heat output in resting muscles. SERCA has been reported to account for ∼25% of the resting metabolism of mouse slow-twitch muscle (Chinet et al., 1992) and 20% of rat fast-twitch muscle (Simonides and Van Hardeveld, 1988), but there are also much lower estimates (e.g., 3–5%; Clausen et al., 1991; Rolfe and Brown, 1997). The rate of heat output from noncontracting human forearm muscles is 0.6 W kg−1 (Table 1), or 600 mW kg−1, so if it is assumed that 25% of this is due to Ca2+ cycling via SERCA, then ATP use by the pump would produce heat at 0.25 × 600 = 150 mW kg−1.
Setting Q to 150 mW kg−1 and solving for gives 0.66 μM. This is probably a lower limit, because would be greater if resting CaM was >50 nM. For the current study, was set to 0.7 μM, which satisfies Eq. 1 and is comparable to estimates for other mammalian fast-twitch muscles (0.3 μM, Vilsen and Andersen, 1992; 0.62 μM, Cantilina et al., 1993).
For SERCA, γ is ∼6 s−1 at 25°C (calculated from data in Lytton et al., 1992; Cantilina et al., 1993; Fortea et al., 2001). was taken to be 711 μM s−1 (Table 2). Setting nSERCA = 2 and rearranging Eq. 2, cSERCA = 711 μM s−1 / (2 × 6 s−1) = 59 μM (expressed relative to myoplasmic water volume). This compares favorably with values for fast-twitch muscles of other mammalian species. Baylor et al. (1983) calculated a pump concentration of 69 μM for rabbit fast-twitch muscle, and the data presented by Ferguson and Franzini-Armstrong (1988) (their Table 2) can be used to calculate a value of 68 μM for guinea pig fast-twitch muscle.
Ca2+-cycling component of resting metabolism constrains the magnitude of SR Ca2+ leak
There is evidence that Ca2+ leaks into the myoplasm from the SR via mechanisms in addition to RYR Ca2+ leak. Ca2+ leaks from the SR, even with the RYRs blocked, and at least part of this leak involves SERCA itself (Simonides and Van Hardeveld, 1988; Macdonald and Stephenson, 2001). A passive, concentration gradient–driven leak from the SR into the myoplasm, with a rate constant was incorporated into the model. The magnitude of Ca2+ leak, by all mechanisms, from the SR into the myoplasm of resting muscle can be estimated from the energy expenditure of SERCA in resting muscle. Assuming that Ca2+ cycling accounts for a heat production of 150 mW kg−1, as calculated above, and that the energy equivalent of ATP turnover is 74 mJ μmol−1 (see above), then SERCA-related ATP turnover would be 150 mJ s−1 kg−1 / 74 mJ μmol−1 ≈ 2.0 μmol kg−1 s−1, and the rate of Ca2+ pumping would be 4 μmol kg−1 s−1. Taking into account of the relative volumes of the myoplasm and SR (Table 3), this uptake rate corresponds to Ca2+ exiting the SR at 69 / 3.8 × 4 = 73 μM s−1 with respect to SR volume. If CaSR at rest is 400 μM (Sztretye et al., 2011b; Manno et al., 2013), this rate of loss is consistent with a rate constant of 73 μM s−1 / 400 μM = 0.18 s−1. This represents the total SR Ca2+ cycling, including leaks via RYRs and any other sources. will, therefore, be 0.18 − s−1, where is the rate constant for Ca2+ efflux from the SR through the RYRs.
Plasma membrane Ca2+ pump (PMCA)
It was assumed that the pump responsible for moving Ca2+ into the t-tubule system is of the PMCA family (Cully et al., 2012). PMCA has a molecular structure similar to SERCA, with the addition of a region enabling kinetic modulation in vivo by calmodulin (Carafoli, 1991), and is localized to the junctional region of skeletal muscle (Sacchetto et al., 1996). The slope of the relationship between Ca2+ pumping rate and [Ca2+] for PMCA is ∼2 (Kosk-Kosicka and Inesi, 1985), indicative of the cooperative binding of 2 Ca2+ in each pump cycle, which is consistent with direct measurements of the pump stoichiometry (Akyempon and Roufogalis, 1982; Valant and Haynes, 1993). It was therefore assumed that the mechanism of PMCA involves the binding and transport of 2 Ca2+ in each ATP-splitting cycle. The sensitivity of PMCA to Ca2+ depends on calmodulin. In the absence of calmodulin, the pump’s Ca50 is 5–10 μM, but this is reduced to <0.5 μM in the presence of calmodulin (Enyedi et al., 1993). EGTA, which was in the solutions used for the fiber studies on which the current analysis is based, also lowers and to about the same extent as calmodulin (Schatzmann, 1973). The maximum pump turnover is similar for SERCA and PMCA. For instance, the maximum rates of ATP splitting by PMCA purified from human erythrocytes (Kosk-Kosicka and Inesi, 1985; Kosk-Kosicka and Bzdega, 1988; Enyedi et al., 1993), measured in the presence of calmodulin, are between 5 and 7 s−1 at 25°C (assuming a Q10 of 2). Therefore, as for SERCA, a maximum turnover of 6 s−1 was used.
Given the structural and biochemical similarities between SERCA and PMCA, the same scheme (Fig. 2 A) was used to model both pumps (Sneyd, 2005; Croisier et al., 2013). Support for this approach can be seen in that reported relationships between the rate of ATP breakdown by PMCA and [Ca2+] can be described well using the model. For example, two empirical relationships between rate of ATP splitting and [Ca2+] are shown in Fig. 3, and both have been fitted with a curve calculated using the PMCA model equations, adjusting only the pump’s Ca2+ affinity. The form of the model predictions matches the data well, further supporting the idea that PMCA behaves in a manner consistent with cooperative binding of 2 Ca2+ in each cycle. The dashed line in Fig. 3 shows the ATPase–Ca2+ relationship produced when the model parameters were optimized to the human muscle fiber PMCA Ca2+ accumulation data of Cully et al. (2018).
Determinants of steady-state CaT: thermodynamics or leak of Ca2+ from t-system?
When the t-system fills with Ca2+, CaT reaches a plateau after ∼60–100 s (Cully et al., 2018). There are two potential causes of this phenomenon. First, the PMCA may reach its thermodynamic limit; that is, the plateau occurs when ΔGCa rises sufficiently, as the concentration gradient increases, to match ΔGATP. An alternative, and not necessarily mutually exclusive, explanation is that there is passive leak of Ca2+ out of the t-system, the rate of which depends on the [Ca2+] gradient between the t-tubule and the myoplasm. Each of these concepts is assessed in the following paragraphs. Predictions of the two models for limiting CaM are compared to experimental data obtained in the presence of tetracaine (Fig. 4, circles); this is the appropriate comparison because in the absence of a Ca2+ leak into the JS, CaJS, which drives the pump, is determined by the Ca2+–EGTA equilibrium; thus, CaJS ≈ CaM and can be estimated with some confidence.
Eq. 5 was used to calculate the relationship between the limit CaT and CaJS for a series of values of ΔGATP; these are shown by the curves in Fig. 4 A. The range of ΔGATP values used for the calculations (40–55 kJ mol−1) span the likely values in the skinned human muscle fibers. In intact skeletal muscle fibers, ΔGATP is typically ∼60 kJ mol−1 (Barclay, 2015), but lower values are likely for skinned fibers, primarily owing to elevated inorganic phosphate concentration (Park-Holohan et al., 2010). The curves in Fig. 4 A can be used to estimate a lower limit for ΔGATP in the skinned fibers. Even if PMCA thermodynamics did not limit Ca2+ accumulation, the pump must be capable of generating the observed Ca2+ gradients. From Fig. 4 A, if ΔGATP was 40 kJ mol−1, then the measured CaT values for CaJS of 28 and 67 nM were approximately fivefold higher than the pump could generate; therefore, ΔGATP cannot have been that low. The curve for ΔGATP of 45 kJ mol−1 is close to the measured data points for CaJS of 28 and 67 nM. This indicates that ΔGATP in the fibers must have been ≥45 kJ mol−1. However, the experimental CaT values for CaJS of 10, 210, and 1,340 nM lie substantially below the limit CaT, so for those values, at least, some other mechanism must limit Ca2+ accumulation in the t-system. In fact, the difference between the forms of the limit CaT–CaJS curves (increasing monotonically) and the experimental data (sigmoidal) precludes a thermodynamic pumping limit determining all the observed CaT values.
To incorporate the t-system Ca2+ leak into the fiber model, it was important to decide the destination of the Ca2+ leaked from the t-system. The leaked Ca2+ could enter either the JS or the myoplasmic space. In trials, it was found that if the Ca2+ leaked into the JS then when CaM was <50 nM, the leak caused a marked increase in CaJS and, consequently, predicted CaT values were higher than those observed. In contrast, when the leak was directed into the myoplasm, there was sufficient buffering to prevent changes in CaM, and the predicted CaT matched the observed CaT–CaM relationship. Therefore, it was assumed that Ca2+ that leaked from the t-system entered the myoplasmic space.
Estimating PMCA concentration
If the rate at which Ca2+ accumulates in the t-system depends on the balance between the rates of Ca2+ uptake into and leak out of the t-system, then three parameters determine t-system Ca2+ accumulation: cPMCA, and Values are not known for any of these parameters. Parameter optimization algorithms can be used to find a combination of values that will allow the model to match experimental data, but in this case such solutions are not unique because cPMCA and are not independent. Therefore, a value for one of the three parameters is required.
cPMCA could be estimated using Eq. 2, with the reported maximum rate of t-system Ca2+ accumulation (using Lagrange’s notation, CaT'′). Cully et al. (2018) reported a mean maximum CaT′ of 170 μM s−1 (with the concentration referring to t-system volume). However, this rate is too high to be consistent with the overall time course of t-system Ca2+ accumulation. For example, when CaM is 67 nM and tetracaine is present, CaT increases with a half-time (t½) between 15 and 20 s, and >60 s is required to reach a final steady CaT (e.g., Fig. 5 A). If the maximum rate of Ca2+ pumping is 170 μM s−1, then by Eq. 2, cPMCA would be 280 μM, with concentration expressed relative to JS volume. Using that concentration, the model predicts that CaT will reach a steady value in ∼15 s, regardless of the value of This suggests that the reported maximum CaT′ is not representative of the overall uptake of Ca2+ into the t-system.
As an alternative, a modeling-based method was used to estimate cPMCA. The simple PMCA/leak model described above (Eqs. 6 and 7) was used to find combinations of cPMCA and that produce records of t-system Ca2+ accumulation consistent with the observed rate of increase in CaT and its final magnitude. Using the best-fit parameters from the preliminary analysis of the pump/leak model (Fig. 4 B) as a starting point, a detailed Monte Carlo analysis was performed in which the time course of increase in CaT was calculated using the model for 5,000 combinations of cPMCA and each selected at random from a range of ±50% of the preliminary values (Fig. 5 B). The criteria set for a record to be considered consistent with observed data were that t½ was between 15 and 20 s and the maximum CaT was between 690 and 770 μM (±5% of the reported mean value; Cully et al., 2018). Records meeting the criteria were identified (i.e., those that passed through the red rectangle and finished between the horizontal red lines in Fig. 5 B), and the combinations of cPMCA and giving rise to such records were plotted (Fig. 5 C). To determine the effect of the third influential parameter, the analysis was performed for three values of that span the range of likely values; is constrained by the range of [Ca2+] corresponding to the sloping section of the CaT–CaM relationship (Fig. 4 B). The results of this analysis are shown in Fig. 5 C.
For each value of ∼200 records were identified as being consistent with the observed time course and magnitude of t-system Ca2+ accumulation; the combinations of cPMCA and that gave those records are plotted in Fig. 5 C. For each value of the range of pump concentrations that can give a time course consistent with those observed is ∼15 μM. As a best estimate of cPMCA, the value in the middle of the complete range was taken; this was 58 μM. This value, with selection of an appropriate is able to accurately simulate the CaT time course across the full range of values tested. cPMCA is expressed relative to the JS volume (Table 3), which is appropriate because the rate of Ca2+ pumping is set by CaJS and the cPMCA in that volume. Note that although cPMCA is numerically similar to cSERCA, when the difference in volume between the JS and myoplasm (Table 3), and in the surface area of the t-system and SR (Eisenberg, 1983), are taken into account, the current analysis suggests that the density of Ca2+ pumps in the t-system membrane is <1% of that for the SR membrane.
In Fig. 5 D, a comparison is shown of the measured record of the time course of CaT from Fig. 5 A and that calculated using the PMCA/leak model with cPMCA of 58 μM, of 0.02 s−1, and of 200 μM (corresponding to pump Ca50 = 40 nM). The time courses match well, although the model indicates that the experimental protocol did not allow sufficient time for CaT to reach a complete plateau. Overall, the analysis presented in this section shows that a model of t-system Ca2+ accumulation consisting of PMCA, driven by a constant [Ca2+], and a concentration gradient-dependent Ca2+ leak from the t-system can accurately describe both the steady-state CaT–CaM relationship measured in the presence of tetracaine (Fig. 4 B) and the time course with which Ca2+ accumulates in the t-system (Fig. 5 D).
Buffering of Ca2+ by EGTA in JS and myoplasm
Buffering of Ca2+ in the SR
In resting muscle fibers, most of the Ca2+ is in the SR, and most of that is bound to Csq and, probably, to other unidentified buffers (Manno et al., 2013). The buffering of Ca2+ in the SR is an important component of the model, as it was assumed that the rate of Ca2+ leak through the RYRs depended on the Ca2+ concentration gradient between the SR and JS. In the protocols used in the fiber experiments, CaJS was maintained in the submicromolar range by its equilibrium with EGTA; in that case, the rate of leak is primarily determined by CaSR (typically 400–1,000 μM) and is thus influenced by the features of Ca2+ buffering in the SR. Additionally, the protocol for measuring RYR Ca2+ leak starts with the SR depleted of Ca2+; therefore, the magnitude of the leak will, at least initially, be time dependent, reflecting the time course with which CaSR increases, and so is likely to be influenced by Ca2+-buffering kinetics.
In Eq. 9, B is buffer, CaB is Ca2+ bound to the buffer molecule, and kon and koff are the rate constants for association and dissociation, respectively. In Eq. 10, BTot is the total concentration of buffer sites, CaB is the concentration of buffer sites with Ca2+ bound, BTot − CaB is the concentration of available buffer sites, and n is the cooperativity coefficient. Manno et al. (2013) estimated that the concentration of Ca2+ buffer sites was ∼70 mM (with respect to SR volume), and this value was used in the model. The rate constants are related to the dissociation constant by Dissociation of the Ca2+-buffer complex is rapid (Baylor and Hollingworth, 2003; Sztretye et al., 2011a), so koff was set to 3,000 s−1. The SR Ca2+ buffering properties of mouse muscle are well described by a sigmoidal relationship, with = 460 μM and n = 3.5 (Manno et al., 2013). With only a slight adjustment to (450 μM), using these values in the current model (solid line, Fig. 6 B) provided a close match to the empirical SR Ca2+ buffering curve described by Manno et al. (2013) (Fig. 6 B, black line). Using these values, = 3,000/e3.5ln(450) = 1.55 × 10−6 μM−1 s−1. Fig. 6 C shows the variation in bound and unbound [Ca2+] as a function of total SR content. A notable feature is that when the SR Ca2+ content is low, the Ca2+ buffering power is low. Consequently, during the initial stages of SR filling (after prior Ca2+ depletion), CaSR will increase with little increase in CaB.
Buffering of Ca2+ by TnC in the myoplasm
TnC binds Ca2+ in its role of linking the contractile response to the release of Ca2+ into the myoplasm. Fast-twitch mammalian muscle myoplasm contains ∼200 μM of TnC Ca2+ binding sites, and at rest, ∼30% of the sites have Ca2+ bound. The question of relevance to the current model is, can TnC can bind sufficient Ca2+ so that the free Ca2+ is not simply that expected from the total Ca2+ in the fiber and the kinetics of Ca2+ binding by EGTA? To address this question, a model was developed of Ca2+ binding by TnC and EGTA. The model is described in detail in Appendix 2. The modeling illustrated that Ca2+ binding by TnC would have a negligible effect on the free Ca2+ calculated from just the Ca2+–EGTA reaction. For example, when the free Ca2+ is intended to be 1 μM, the presence of TnC would reduce the actual free [Ca2+] by 0.1%; the effect is even smaller for lower free [Ca2+]. Consequently, it was not necessary to include Ca2+ binding by TnC in the model.
Diffusive exchange of Ca2+ between JS and myoplasm
In the absence of a leak of Ca2+ into the JS, there is little difference between CaJS and CaM because both are controlled by the equilibrium with EGTA. In the presence of a leak, it would be expected that CaJS increases as Ca2+ is added to the JS, despite the buffering by EGTA. The observation that a steady, submaximal CaT is achieved when CaM is at levels too low to achieve maximal pumping by PMCA indicates that there is a mechanism that limits the increase in CaJS. We assumed that this was diffusion of Ca2+ between the JS and myoplasm. The diffusion rate was adjusted so that the model steady-state CaT matched the experimentally determined values. That is, in the steady state, the sum of the leak of Ca2+ from the t-system into the myoplasm and the diffusive movement from the JS to the myoplasm equaled the rate at which Ca2+ was added to the JS via the RYRs. The diffusive exchange was modeled as a simple movement of Ca2+ depending on the difference (CaJS − CaM).
Relative volumes of physical compartments
The physiologically relevant Ca2+ concentrations (e.g., for stimulating pumps or driving diffusive exchanges) are those expressed with respect to the volume of water within individual compartments rather than to fiber water volume. Therefore, the relative volumes of the different compartments were taken into account when determining the concentration changes resulting from exchanges between compartments. For example, in human vastus lateralis muscle, the relative volume of the SR is ∼100 times that of the JS (Eisenberg, 1983), so that movement of Ca2+ from the SR to the JS results in an increase in [Ca2+]JS that is 100-fold greater than the corresponding decrease in [Ca2+]SR. The proportions of fiber volume accounted for by each physical compartment are shown in Table 3.
The JS is a conceptual rather than physical space, and its volume is difficult to both define and estimate; it is a region between the RYRs and the adjacent t-tubule. We estimated VJS, relative to that of the t-tubule (VTT), on the basis that (1) t-tubule diameter is 85 nm (Jayasinghe and Launikonis, 2013), (2) the distance between the SR and t-tubules is 12 nm (Mitchell et al., 1983; Franzini-Armstrong et al., 1999), (3) the RYRs are 30 nm in diameter and the spacing between the centers of closest adjacent units is 100 nm (Franzini-Armstrong et al., 1999), and (4) the RYRs are arranged in a double layer, with one layer on either side of the t-tubule (Jayasinghe et al., 2014). The volume of a space with the same diameter as the RYR and the t-tubule was calculated, and it corresponded to ∼4% VTT if the t-tubule cross-section is circular and 5% if elliptical (Jayasinghe and Launikonis, 2013), with the major axis length twice that of the minor axis. This range of values is consistent with the 4.5% calculated by Dulhunty et al. (1992). In the simulations, VJS was set to 5% VTT.
Unless otherwise specified, concentrations are expressed relative to the compartment in which the relevant molecular species is contained. For example, concentrations in the myoplasm are expressed with respect to myoplasmic water volume. The data used to calculate myoplasmic water volume and the concentration conversion values are presented in Table 3.
Composition of myoplasmic and JSs
The rates of Ca2+ movement by SERCA and PMCA depend on the concentrations of ATP and its metabolites and on pH. These values are largely determined by the composition of the solutions used with the skinned fiber preparation (Cully et al., 2018). Assumed concentrations (in mM) were as follows: myoplasm/JS: 8 ATP, 0.05 ADP, and 1 Pi; t-tubule: Ca2+, time dependent and calculated using the model. pH in all spaces was taken to be 7.1. The rate of Ca2+ pumping by SERCA and PMCA is strongly influenced by ADP concentration and pH (Inesi and Hill, 1983; Tran et al., 2009). It was assumed that [ADP] and pH were constant during an experiment, with pH buffered by HEPES in the bathing solution and ATP, and hence ADP, buffered by endogenous creatine kinase, presumably bound to intracellular structures (Dzeja and Terzic, 2003).
The experimental data used in this study are those referred to in Cully et al. (2018) that describe the uptake of Ca2+ into the t-system of fibers from human vastus lateralis muscle. The data are from the control group in that study and as such had no apparent abnormalities in Ca2+ handling. Briefly, the experimental protocol involved first depleting the SR and t-system of Ca2+ by exposing the fibers to caffeine, which kept the RYRs and store-operated Ca2+ channels on the t-system open. The caffeine was removed, and the fibers were placed in a solution with a known Ca2+ concentration. The time course with which Ca2+ entered the t-system was monitored using the Ca2+-sensitive compound rhod-5N trapped inside the t-tubules.
Implementing the model
The modeling was carried out using Maple (v2016.1) and MapleSim (v2016) software (Maplesoft). The model equations were deployed, and parameter values were entered in the general mathematics environment provided by Maple. The system of differential equations was solved numerically using the Rosenbrock method (Rosenbrock, 1963). To be consistent with the experimental protocol (Cully et al., 2018), simulations commenced with the SR and T systems depleted of Ca2+ so that all the Ca2+ was in the myoplasm and JS and in equilibrium with EGTA. The total Ca2+ was set to give the desired myoplasmic free [Ca2+], as determined by the EGTA–Ca2+ equilibrium (Fig. 6 A). The differential equations are shown in Appendix 1. Once the model was developed, the equations were imported into MapleSim to create a fast, compiled version of the model, which was used for subsequent analysis using Maple. Using a computer with a 2.2-GHz processor and 8 GB of memory, the compiled version simulated 120 s of t-system Ca2+ accumulation in 8 ms.
To quantify the rate of Ca2+ leak from RYRs, cPMCA was fixed at 58 μM, and and were treated as variables whose values were those that provided the least-squares fit of the model to the empirical CaT–CaM relationship determined in the presence of tetracaine (see Results). This was done by creating a function that returned the model-calculated maximum CaT for specified values of CaM, and and using that function in a nonlinear regression fitting routine implemented using Maple software.
Identifying parameters that influence t-system Ca2+ accumulation
Parameters that influenced t-system Ca2+ accumulation were identified using a Monte Carlo analysis in which the model was solved for many random combinations of parameter values. The analysis was performed using nine parameters selected on the basis that their location in the model might allow them to influence t-system Ca2+ accumulation. The results of the analysis are summarized in Table 4 and Fig. 7. In Fig. 7, influences of four parameters are illustrated graphically: the cloud of light gray symbols are the CaT (Fig. 7, A–D) and t½ (Fig. 7, E–H) values calculated from 4,000 CaT versus time curves with random combinations of the nine parameters, and the solid lines show the variation in CaT or t½ as a function of only the parameter of interest. To quantify the influence of each parameter, the square of the correlation coefficient (r2) was calculated from the least-squares fit of a second-order polynomial through the data cloud, with the parameter of interest as the independent variable. r2 is the fraction of the total variance, arising from all the parameters assessed, that can be attributed to the parameter of interest. The r2 values for all the parameters assessed are shown in Table 4, which includes analyses for CaM of both 67 and 1,340 nM.
The four parameters with the greatest influence on Ca2+ accumulation in the t-system were those directly involved in the t-system—PMCA concentration and Ca2+ sensitivity (quantified by and the magnitude of the Ca2+ leak from the t-system—and the Ca2+ leak through the RYRs. These parameters accounted for ∼75% of the total variance in steady-state CaT and t½. The extent of t-system filling was mainly determined by cPMCA (Fig. 7 A) and t-system leak rate (Fig. 7 C); this is consistent with the notion that the level of filling reflects the balance between rate of Ca2+ uptake into the t-system, which depends on cPMCA, and the leak of Ca2+ from the t-system. The size of the RYR leak rate constant also influences the magnitude of CaT, which is desirable given that the purpose of the protocol was to use changes in CaT to estimate RYR leak. It should be noted, however, that although influenced the extent of Ca2+ accumulation when CaM is lower than that required to fully activate PMCA (e.g., 67 nM), it had no influence when CaM is high (1,340 nM). In the latter case, the pump is operating maximally even in the absence of RYR leak, so any further addition of Ca2+ to the JS by the leak would not affect the rate at which Ca2+ is pumped into the t-system.
It is apparent from Table 4 that parameters associated with the SR Ca2+ pump (the right three columns of Table 4) had little influence on Ca2+ accumulation in the t-system. This indicates that the time course of increase in CaSR, and any effect that has on RYR Ca2+ leak and thus CaJS, does not limit the pumping of Ca2+ from the JS into the t-system.
Steady-state CaT–CaM relationship
Uptake of Ca2+ into the t-system has been characterized by the relationship between the ultimate extent of Ca2+ accumulation and the concentration of Ca2+ in the myoplasmic solution (i.e., the CaT–CaM relationship; Fig. 4 B). From Fig. 7, it is apparent that steady-state CaT depends on cPMCA, and the rate constants for the Ca2+ leaks from the t-system and through the RYRs. The value of cPMCA was fixed (see Materials and methods) but the values of and are not known. Those values were determined using the model by finding the parameter values that gave the least-squares fit of the model to the experimental CaT–CaM data. There were two steps to this analysis. First, and were determined using data obtained in the presence of tetracaine, because in that case it was assumed there was no RYR Ca2+ flux, so = 0. It was then assumed that was unaffected by tetracaine but that, as preliminary analysis revealed, did depend on whether the RYRs were functional. Therefore, the second part of the analysis was to use data obtained without tetracaine to determine (in the absence of tetracaine) and The results of the fitting procedure are shown in Fig. 8 A, and the parameter values are given in Table 5.
The shape of the CaT–CaM relationship predicted by the model matches the experimental data well (Fig. 8 A), indicating that the model included sufficient detail to describe t-system Ca2+ accumulation in the experimental preparation. In the absence of RYR Ca2+ leak, the calculated and values (Table 5) give rise to a CaT–CaM relationship with a of ∼67 nM. Withdrawal of tetracaine had two clear effects. (1) The leak of Ca2+ from the t-system was reduced, which increased steady-state CaT values. in the absence of tetracaine was 60% of that in the presence of tetracaine. (2) Ca2+ leak through the RYRs was apparent = 0.035 s−1), which shifted the CaT–CaM curve to the left of that determined with tetracaine; was reduced to 42 nM. It is difficult to distinguish these two effects—increased RYR leak and decreased t-system leak—by inspection, but they were readily revealed and quantified using the model curve-fitting procedure.
The rate constants for the two Ca2+ leaks (Table 5) were combined with known values for CaSR and CaT to calculate the rates of Ca2+ exchanges among fiber compartments. These calculations were done using steady-state data corresponding to CaM of 67 nM, this being approximately the free [Ca2+] in unstimulated, intact muscle fibers, and assuming that RYR leak was as determined in the skinned fiber preparation in the absence of tetracaine. The steady-state movements of Ca2+ in the skinned fiber preparation are summarized in Fig. 9. The flows were calculated using concentrations of Ca2+ in the SR and t-system and the rate constants for leaks from those compartments (Table 5). For example, the rate of Ca2+ leak from the t-system can be calculated: × steady-state CaT = 0.013 s−1 × 1,200 μM = 15.6 μM s−1, and this is the same as the rate of uptake of Ca2+ into the t-system by PMCA in the steady state. This sets the rate at which Ca2+ is removed from the JS, and so on. In Fig. 9, the flows are quantified according to the volumes of both the origin and destination compartments (Table 3). For example, an outflow (indicated by −) of 14 μM s−1 from the SR produces an influx (+) into the JS, in terms of JS volume, of 14 × (3.8 / 0.014) = 3,800 μM s−1. Note that in a steady state, the amount of Ca2+ bound to EGTA will be constant and thus does not need to be considered in this analysis.
A notable feature of the Ca2+ flows is that only a small fraction of the Ca2+ calculated to leak into the JS is pumped into the t-system: of the influx of 3,800 μM s−1, 312 μM s−1, or 8%, is pumped into the t-system, and the remainder must diffuse into the myoplasmic space. This is, at least in part, a consequence of the low density of PMCA in the JS, as calculated in Materials and methods. A second notable aspect is that ∼20% of the Ca2+ pumped into the SR leaks into the JS, and the remainder leaks into the myoplasm. Therefore, under resting conditions in the skinned fiber, most of the cycling Ca2+ leaked from the SR into the myoplasm, presumably via a SERCA-related mechanism; one-fifth leaked into the JS via RYRs; and only a small fraction of that passed through the PMCA.
Time course of Ca2+ movements during t-system Ca2+ accumulation
The model provides a means of estimating the time course of Ca2+ movements during the protocol used to determine Ca2+ t-system accumulation. In Fig. 8 B, the time courses of Ca2+ concentration in the myoplasm and JS are shown for a RYR leak of the magnitude given in Table 5. In the absence of a leak of Ca2+ into the JS through the RYRs, CaJS stayed close to its initial value of 67 nM throughout the period over which the t-system accumulated Ca2+ (blue line, Fig. 8 B). That constancy reflected the buffering of Ca2+ by EGTA and diffusion of Ca2+ into the JS from the myoplasm. In the presence of a leak of Ca2+ into the JS, CaJS increased toward ∼100 nM (red line, Fig. 8 B). In Fig. 8 C, the time courses are shown of Ca2+ entering the SR and t-system and binding to buffers in the SR. Note that when CaM is 67 nM, the rate at which SERCA operates is ∼1% of its maximum, so the filling of the SR is relatively slow. In the SR, most Ca2+ binds rapidly to the SR buffers (magenta curve, Fig. 8 C). The buffer characteristics allow CaSR to increase rapidly at the start of filling (green curve, Fig. 8 C), thus stimulating the RYR Ca2+ leak when tetracaine is not present. The t-system begins to fill immediately, fueled initially by the free Ca2+ in equilibrium with EGTA in the JS when the protocol begins. In the absence of tetracaine, t-system filling is enhanced by the rise in CaJS due to the RYR leak (compare with tetracaine [blue curve] to without tetracaine [red curve]; Fig. 8 C).
How precisely can the RYR Ca2+ leak be quantified?
A benefit of using the curve-fitting approach to quantifying the RYR leak is that the parameters derived from the fitting procedure can be used to calculate confidence intervals (CIs) for the model parameters whose values were determined from the fitting procedure. Taking the data obtained in the absence of tetracaine as an example and using the raw data points when fitting the model to the data (Fig. 10 A), a confidence contour can be calculated (Motulsky and Christopoulos, 2004). The confidence contour is an ellipse (outlined by open circles, Fig. 10 B) that encompasses combinations of parameter values that would give rise to a curve that would not differ significantly, at the 5% level in this case, from the best-fit curve. The parameter combination that gave the best fit is indicated by the solid symbol in the middle of the ellipse (Fig. 10 B, filled circle). The extent of the excursion of the ellipse gives the 95% CI for each parameter, indicated by the dashed lines.
The 95% CI for the RYR leak rate constant is 0.007–0.041 s−1. The interval excludes 0, which indicates that the method used to quantify the RYR leak was sufficiently sensitive to detect the leak that occurred in the human skinned fibers. The 95% CI would likely be narrowed if data points were obtained at additional CaM values on the sloping part of the CaT–CaM curve.
The 95% CIs for can be determined in the same way (Fig. 10 B and Table 4). In this case, the 95% CIs for the data with and without tetracaine do not overlap, so it can be concluded that the t-system Ca2+ leak rate differed significantly depending on whether the RYRs were functional. That is, when Ca2+ entered the JS through the RYRs, the leak of Ca2+ from the t-system was reduced compared with that occurring when CaJS was held constant by the EGTA buffering.
Plotting CaT versus CaM or CaJS?
The convention adopted by Launikonis and colleagues (Cully et al., 2016; Cully et al., 2018) is to plot the extent of t-system Ca2+ accumulation as a function of CaM. The t-system filling is driven by the PMCA and so must reflect CaJS rather than CaM. However, CaJS is unknown; hence the practice of using CaM as the independent variable. The model provides the means to estimate CaJS, and in Fig. 11 A, a comparison is shown of steady-state CaT values plotted against CaM and the model-derived CaJS. For this comparison, the effect of the change in has been removed from the comparison by normalizing steady-state CaT values obtained with and without tetracaine data by their respective maxima (i.e., values at CaM of 1,340 nM). When plotted as a function of CaM, the effect of withdrawing tetracaine, thus allowing the RYR Ca2+ leak, is to shift the normalized CaT–CaM relationship to the left. This effect is eliminated when CaT is plotted against CaJS; in that case, the two curves (solid blue and dashed red lines, Fig. 11 A) superimpose. Remembering that changes in have been excluded from this analysis, the observation that the CaT–CaJS curve overlies the CaT–CaM curve in the presence of tetracaine is important because it confirms (1) that the PMCA properties are unaffected by the presence of tetracaine, (2) that t-system filling simply reflects the combined properties of PMCA and the t-tubule Ca2+ leak, and (3) that the shift of the CaT–CaM curve to the left when tetracaine is removed from the fibers reflects the increase in CaJS in the presence of RYR Ca2+ leak. The last point highlights the benefit of using CaM as the independent variable: the size of the horizontal difference between the normalized curves reflects the magnitude of the Ca2+ leak into the JS, as illustrated in Fig. 11 B.
Applying the model to data from fibers with malignant hyperthermia (MH)-causative RYR mutation
Cully et al. (2018) measured CaT–CaM curves for human muscle fibers exhibiting a MH-causative RYR mutation. The model was used to quantify RYR Ca2+ leak for these fibers. For this analysis, it was recognized that the steady-state CaT value measured without tetracaine at the highest CaM (indicated by * in Fig. 12) is much lower than the point at the next lowest CaM. Cully et al. (2018) suggested that this reflected activation of store-operated calcium entry (SOCE), which occurred only at the highest CaM used. No SOCE mechanism was included in the model, so the unusual data point was excluded from the curve-fitting procedure. The model analysis indicated that was twofold higher in the MH fibers (385 μM, corresponding to = 72 nM) than in the control fibers (Table 5) and that was 0.08 s−1, more than twice as great as that for fibers without the MH pathology (Table 5). If SR function is otherwise unaltered by the RYR mutation, then the model indicates that CaSR would be little affected by the greater RYR leak (i.e., SR Ca2+ is well buffered). The absolute Ca2+ leak, assuming CaSR is 400 μM, would be ∼32 μM s−1. Thus, RYR Ca2+ fluxes in resting muscle fibers are not fixed but rather can be modulated by functional alterations to the RYRs.
The model accurately predicted Ca2+ accumulation by the t-system in the skinned fiber preparation. It predicted the form of the CaT–CaM relationship (Figs. 4 B and 8 A), the effect of tetracaine on that relationship (Fig. 8 A), and the time course with which Ca2+ accumulates in the t-system (Fig. 5 D). This was achieved with most parameter values determined from information in the literature. The only parameters that required fine-tuning to align the model to the experimental data were the rate constants for Ca2+ leaks from the t-system into the myoplasm and from the SR into the JS and the Ca2+ sensitivity of the PMCA, for which values were unknown. The experimental design used by Cully et al. (2018) sought to identify the presence of a Ca2+ leak through the SR Ca2+ release channels by comparing the steady-state CaT–CaM relationships obtained in the presence and absence of the RYR inhibitor tetracaine. From simple observation of the CaT–CaM curves with and without tetracaine, it was difficult to clearly establish whether the relationship altered as would be expected if there was a flow of Ca2+ through the RYRs, because a decrease in the t-system Ca2+ leak alone would also shift the relationship to the left. However, with the available experimental data, the statistical sensitivity of the model-based determination of the RYR leak was adequate to establish that there was a leak through the RYRs in the absence of the RYR blocker, and its magnitude was calculated. Furthermore, use of the model enabled the rate of Ca2+ leak from the t-system to be quantified, in both the presence and absence of RYR leak.
It was important to ensure that the calculated magnitude of the leak was not dependent on the choice of parameter values. The sensitivity analysis indicated that the parameters with most influence on t-system Ca2+ accumulation were, perhaps not surprisingly, those directly associated with the t-system. The value of only one of those, PMCA concentration, was set at a predetermined value. If the pump concentration was different from that assumed, fitting the model to the data would give different values for and but the value of would be unaltered.
JS volume influences the estimated value of as indicated by the sensitivity analysis (Table 4). JS volume is difficult to quantify because the region is a conceptual, rather than physically defined, space. Its assumed volume is important in the model in that it determines how much CaJS increases for a given RYR Ca2+ leak. If a larger JS volume is assumed, then a faster RYR Ca2+ leak is required to produce the observed steady-state CaT. The magnitude of the effect, however, is not great: if JS volume was 50% greater than the value used, then determined from the curve fitting would be increased by ∼12%. Overall, our estimate of seems fairly robust.
A framework for interpretation of experimental data
The analysis presented in Results provides a framework for interpreting the CaT–CaM relationship. In Fig. 13, the effects of a RYR leak, a change in t-system Ca2+ leak and the combined effects are illustrated. The presence of a RYR Ca2+ leak is evident through a pure leftward shift of the relationship determined in the absence of tetracaine compared with that in the presence of tetracaine (Fig. 13 A). Changes in the rate of leak of Ca2+ from the t-system upon the removal of tetracaine are evident as a vertical shift in steady-state CaT values (Fig. 13 B). If decreases, steady-state CaT values increase, and conversely, if increases, CaT values will decrease. When there is both an RYR leak and a change in t-system Ca2+ leak, the two effects can be difficult to distinguish by observation alone (Fig. 13 C). The model curve-fitting procedure described here allows the RYR and t-system Ca2+ leaks to be identified and quantified.
Measurements of t-system filling with CaM set to 1,340 nM provide important information. That Ca2+ concentration exceeds that required to fully activate PMCA. As shown in Fig. 3 B (curve c), the rate of Ca2+ pumping by PMCA reaches its maximum at ∼400 nM. Therefore, if CaM is >400 nM, any increase in CaJS due to RYR leak will have no effect on the rate of Ca2+ pumping by PMCA or on steady-state CaT. This is evident from the absence of correlation between and steady-state CaT when CaM is 1,340 (Table 4). Thus, the tetracaine dependence of steady-state CaT measured at 1,340 nM reflects changes in only unsullied by an increase in CaT shows that has decreased, whereas a decrease shows that has increased. This is important because, as described earlier (see Fig. 12), Cully et al. (2018) showed that in fibers from people with a susceptibility to MH, steady-state CaT at CaM of 1,340 nM was markedly reduced when measured in the absence of tetracaine compared with that when measured in the presence of tetracaine. This is consistent with a substantial increase in most likely via activation of SOCE, independent of any changes in RYR leak, at least when CaJS is high.
How much could RYR Ca2+ leak contribute to resting metabolism?
The cycling of Ca2+ between the SR and myoplasm requires SERCA to move the Ca2+ into the SR. SERCA consumes ATP, which will be regenerated by oxidative phosphorylation in the mitochondria. In a resting muscle in a metabolic steady state, the net heat production arises from the oxidative regeneration of ATP used by SERCA. The SERCA-related heat production associated with the calculated Ca2+ fluxes shown in Fig. 9 can be calculated on the basis of a Ca2+:ATP stoichiometry of 2 and the production of 74 mJ (μmol ATP)−1 (see Introduction). If resting CaM is stable (Carroll et al., 1995) in the face of a continual flux of Ca2+ from the SR into the myoplasm, then Ca2+ must be removed from the myoplasm at the same rate that it enters. From Fig. 9, the RYR efflux of Ca2+ from the SR was 14 μM s−1. To balance that flow, taking account of the relative volumes, this would require SERCA to remove Ca2+ from the myoplasm at 0.77 μM s−1, using ATP at half that rate, 0.39 μM s−1. The associated rate of heat production would be 74 mJ μmol−1 × 0.39 μmol liter−1 s−1 = 29 mW kg−1. With the scheme shown in Fig. 9, 8% of the RYR Ca2+ leak is pumped into the t-system before leaking into the myoplasm and being returned to the SR via SERCA; that is, this component is pumped twice in each Ca2+ leak/uptake cycle. This would increase the heat production associated with the RYR Ca2+ leak by 8%, to 31 mW kg−1.
In comparison, the estimated resting heat production in human forearm muscle is 600 mW kg−1 (Table 1), so the RYR Ca2+ leak component is 100 × 31/600 = 5% of muscle resting heat production. The magnitude of this value is consistent with an experimental determination. Chinet et al. (1992) determined the contribution of PMCA to resting heat production of mouse slow-twitch muscle by measuring the decrease in muscle heat production when PMCA was blocked by pharmacological inhibition of calmodulin. It can be difficult to interpret the results of experiments involving blocking part of the Ca2+ cycling process, but their results were similar to the current analysis: PMCA accounted for 3.2% of resting heat production. Together with the current estimate, this supports the idea that under normal physiological conditions and in the absence of pathological alterations in Ca2+ handling, RYR Ca2+ leak makes only a small contribution to skeletal muscle resting metabolism. However, it is clear that RYR Ca2+ flux can vary across a huge range. The Ca2+ leak calculated in Results for MH fibers was 2.3 times greater than that for fibers with no RYR pathology and thus could be expected to produce ≥60 mW kg−1, which could increase muscle resting metabolism by 5%. At the other extreme to the leak in resting fibers, in response to stimulation, SR [Ca2+] decreases at rates of up to 50 mM s−1 (Sztretye et al., 2011a), >3,000 times greater than the leak measured in resting fibers.
Non-RYR Ca2+ leak from the SR
The experimental evidence presented by Chinet et al. (1992) indicated that ∼25% of resting metabolism is related to SERCA activity. The simplest interpretation of this observation is that there is a continual leak of Ca 2+ from the SR into the myoplasm or JS that is pumped back into the SR with the metabolic cost (as either heat production or O2 consumption) reflecting the mitochondrial activity required to support SERCA-related ATP turnover. The current analysis shows that RYR Ca2+ leak is only one-quarter of the total leak required to account for the SERCA-dependent Ca2+ leak (Fig. 9). If this is correct, what accounts for the majority of the leak? Ca2+ leaks from the SR of resting fibers (Macdonald and Stephenson, 2001; Lamboley et al., 2014) and from SR vesicles (e.g., Simonides and Van Hardeveld, 1988) via a mechanism directly involving SERCA. These leaks were not associated with the normal pump cycle and thus are not related to pump reversal or uncoupling (both unlikely under physiological conditions of low myoplasmic Ca2+ and SR free [Ca2+] of ∼500 μM). In fact, a leak mechanism, independent of the ATP-splitting pump cycle, seems to be a characteristic of SERCA (Jilka et al., 1975); Berman (2001) has described this as “channel-like behavior.”
Lamboley et al. (2014) reported that the SERCA-related Ca2+ leak in skinned human fibers was 6–8 μM s−1 (recalculated from the original data to express concentrations with respect to myoplasmic volume, allowing direct comparison with values in Fig. 9). That was measured at pH 8.5, which those authors suggested would increase the leak compared with that at physiological pH. In that case, those values are probably in reasonable agreement with the 3.2 μM s−1 from the current analysis (Fig. 9) that arises from the assumption that 25% of resting metabolism arises from Ca2+ cycling. Therefore, it would appear that assumption is not unreasonable, and given the relative magnitudes of the RYR and SERCA Ca 2+ leaks, a comprehensive investigation of the role of Ca2+ cycling in resting metabolism of skeletal muscle should include both leak pathways.
The model developed for this study can accurately predict the characteristics of t-system Ca2+ accumulation and can be used to quantify the Ca2+ leak through the RYRs. The main analysis used data from human muscle fibers with no Ca2+ handling–related pathology. In these fibers, there is only a small flux of Ca2+ through the RYRs when the fibers are unstimulated. However, RYR Ca2+ fluxes can potentially vary greatly. For example, in resting fibers, larger RYR fluxes occur in pathological states, such as MH (Cully et al., 2018), and in muscle from mice developed to study variations in RYR function (Lamboley et al., 2021). It is apparent from the CaT–CaM relationships published in those papers that there are changes with genotype in not only but also and The current model provides a means by which these changes can be distinguished and quantified.
The heat generated, mostly by SERCA, in balancing the small RYR leak is only a small fraction of the resting metabolic rate of skeletal muscle. Although Ca2+ pumped into the t-system by PMCA is used by Launikonis and colleagues as a biosensor, it is unlikely that PMCA contributes greatly to muscle heat production. As illustrated in Fig. 9, because of the relatively small amount of PMCA compared with SERCA, the majority of Ca2+ that leaks into the JS will be removed from the myoplasm into the SR via SERCA, and that will be the source of most of the heat production. In the current context, the discrete localization of PMCA in the JS, allied with the protocol of comparing t-system Ca2+ accumulation in the presence and absence of an RYR inhibitor and a quantitative model, makes for an extremely sensitive sensor of Ca2+ movement through the RYR in resting muscle.
The current model was a development of a previous one that was used to predict Ca2+ movements and heat production in intact muscle fibers in response to stimulation (Bakker et al., 2017; Barclay and Launikonis, 2021). The improvements in the current model include an enhanced spatial representation of the fiber, a more physiological model of SERCA, and the inclusion of PMCA. There are more sophisticated models of cellular Ca2+ handling (e.g., Baylor and Hollingworth, 2007), but a strength of the current approach is the extension from the microscopic Ca2+ movements within a cell to the macroscopic consequence of those events in terms of heat production. Incorporation of heat production in the model is useful both for comparison with observed values and as a constraint on the parameters of SERCA. Some models of Ca2+ movements incorporate diffusion of mobile Ca2+ entities within the fiber (Cannell and Allen, 1984; Baylor and Hollingworth, 1998). This was deemed unnecessary in the current model because the measured variable, the time course of t-system filling, occurs on a scale of tens of seconds, whereas diffusion on an intrafiber scale (i.e., over micrometer distances) occurs on a millisecond time scale. Any effects of diffusion would be obscured by the slow removal of Ca2+ from the JS by PMCA.
Ultimately, it is the function of intact muscle fibers that is of interest in relation to the heat produced by resting and contracting muscle fibers. Integrating the current skinned fiber model of Ca2+ movements in resting fibers with our previous model of Ca2+ movements in intact fibers in response to stimulus-induced Ca2+ release would be a valuable step on the path to integrate data from skinned fibers into a model that can relate intracellular Ca2+ movements to the metabolism of resting and active intact muscle fibers.
Model equations and parameter values
The skinned fiber model depicted in Fig. 1 was described mathematically as a set of ordinary differential equations. The equations are listed in Table A1. They are categorized as buffering reactions, those involving EGTA or buffers in the SR, leaks, and combined equations. The latter incorporate all movements in or out of a physical compartment and include corrections for the change in volume moving from one compartment to another. Details of the model used for Ca2+ pumping by SERCA and PMCA are given in the following section.
Model parameter values
The model used for SERCA and PMCA was a two-state model (Fig. 2 A) developed by Tran et al. (2009). That was based on an eight-state biochemical model (Fig. A1). The apparent rate constants in the two-state model are functions of the rate constants and dissociation constants in the eight-state scheme as indicated in Fig. A1. In normal operation, the reactions proceed in the clockwise direction with fully cooperative binding of two Ca2+ from the cytoplasm (from state P2 to P4) that are transported across the membrane (P5 to P6) and released into the SR (P6 to P8). Each cycle involves the binding and hydrolysis of ATP, with products ADP and Pi being released at the steps indicated. The model involves H+ both as part of Ca2+/H+ countertransport and in competition with Ca2+ binding. The parameter values for the Ca2+ pump used in the current model are mostly those derived by Tran et al. (2009), with the exception of the Ca2+ dissociation constant on the myoplasmic side of the pump (Kd,Cai in Fig. A1), which was adjusted to match the empirical Ca50 values for SERCA and PMCA (see Materials and methods), and which was adjusted to match the maximum turnover of the pump to that determined empirically (see Materials and methods). Values of all the parameters for the Ca2+ pumps, buffer reactions, and concentrations of cellular constituents are given in Table A2.
Ca2+ binding by EGTA and TnC
An analysis was carried out to see whether Ca2+ binding by TnC would affect Ca2+–EGTA buffering to the extent that free Ca2+ calculated from considering Ca2+–EGTA without taking account of TnC will be inaccurate.
In Fig. A2 A, the calculated time course of the (theoretical) increase of CaM as Ca2+ dissociates from EGTA is shown. The time taken for this to occur is affected by TnC, being slower (on a millisecond scale) when TnC is present. The important result is that the final concentration of unbound Ca2+ was not visibly altered by the inclusion of TnC in the model. By comparing the time courses of CaM, CaTn, and Ca2Tn, it can be seen that EGTA effectively buffers CaM even while Ca2+ is binding to TnC (Fig. A2 B). The reason for this is that the Ca2+–EGTA reaction is faster than Ca2+ binding by TnC, so free Ca2+ that binds to TnC is rapidly replaced by Ca2+ dissociating from EGTA.
Close inspection of calculated steady-state CaM values reveals very small differences in CaM depending on whether the effects of TnC are considered (Table A5). CaM is lower when TnC is considered than when it is ignored. However, the magnitude of the difference is small for the range of CaM values used in the current study. For example, when the anticipated CaM is 67 nM, the actual CaM, taking account of the effects of TnC, would be 0.006% lower. The effect is larger for higher values of CaM, reaching 0.1% when CaM is 1 μM. For practical purposes, these effects are negligible.
Eduardo Ríos served as editor.
This work was supported by Australian Research Council Discovery Projects grants DP180100937 and DP200100435 and AFM Telethon Research Grant 22181 (to B.S. Launikonis).
The authors declare no competing financial interests.
Author contributions: C. Barclay: Conceptualization, formal analysis, methodology, software, writing—original draft, review and editing. B. Launikonis: Conceptualization, data curation, funding acquisition, project administration, writing—review and editing.
This work is part of a special issue on excitation–contraction coupling.