Comprehensive Simulation of Ca2+ Transients in the Continuum of Mouse Skeletal Muscle Fiber Types

Mag-Fluo-4 has revealed differences in the kinetics of the Ca2+ transients of mammalian fiber types (I, IIA, IIX, and IIB). We simulated the changes in [Ca2+] through the sarcomere of these four fiber types, considering classical (troponin –Tn–, parvalbumin –Pv–, adenosine triphosphate –ATP–, sarcoplasmic reticulum Ca2+ pump –SERCA–, and dye) and new (mitochondria –MITO–, Na+/Ca2+ exchanger –NCX–, and store-operated calcium entry –SOCE–) Ca2+ binding sites, during single and tetanic stimulation. We found that during a single twitch, the sarcoplasmic peak [Ca2+] for fibers type IIB and IIX was around 16 µM, and for fibers type I and IIA reached 10–13 µM. The release rate in fibers type I, IIA, IIX, and IIB was 64.8, 153.6, 238.8, and 244.5 µM ms−1, respectively. Both the pattern of change and the peak concentrations of the Ca2+-bound species in the sarcoplasm (Tn, PV, ATP, and dye), the sarcolemma (NCX, SOCE), and the SR (SERCA) showed the order IIB ≥ IIX > IIA > I. The capacity of the NCX was 2.5, 1.3, 0.9, and 0.8% of the capacity of SERCA, for fibers type I, IIA, IIX, and IIB, respectively. MITO peak [Ca2+] ranged from 0.93 to 0.23 µM, in fibers type I and IIB, respectively, while intermediate values were obtained in fibers IIA and IIX. The latter numbers doubled during tetanic stimulation. In conclusion, we presented a comprehensive mathematical model of the excitation–contraction coupling that integrated most classical and novel Ca2+ handling mechanisms, overcoming the limitations of the fast- vs. slow-fibers dichotomy and the use of slow dyes.


Introduction
In mammalian skeletal muscle fibers, the action potentials (AP) lead to contractions mediated by the release of Ca 2+ from the sarcoplasmic reticulum (SR), in a process known as excitation-contraction coupling (ECC) [1]. After release, the Ca 2+ ions bind to a diversity of sites, which include troponin (Tn), parvalbumin (PV), and adenosine triphosphate (ATP). They also flow into the mitochondria (MITO) [1], before being transported back to the SR by a Ca 2+ ATPase (SERCA). Small amounts of Ca 2+ are transported outside the cell by the Na + /Ca 2+ exchanger (NCX). Store-operated Ca 2+ entry (SOCE) allows Ca 2+ enter the fiber through Orai1, as a response to the intra SR Ca 2+ sensing function of the stromal interaction molecules (STIM) [2][3][4]. In fast twitch fibers, SOCE acts in a transient fast mode during an individual AP and after each AP in a train of stimulations [5].
Experiments in muscle fibers loaded with fluorescent indicators have revealed sizeable differences in the kinetic parameters of electrically elicited Ca 2+ transients according to the continuum of fiber types [6]. These variances arise due to the differential quantity and kinetics of the proteins and organelles involved in Ca 2+ handling present in the muscle fibers [1,6]. For instance, SERCA is more abundant in fast fibers [7,8], and PV is almost negligible in slow, but abundantly found in the fastest fibers [9].
Mathematical models that integrate information obtained on mammalian ECC under different experimental conditions have been used to simulate Ca 2+ transients in both slow-twitch and fast-twitch fibers [10,11]. However, the nature of the fibers used has not always been doubtlessly established. Furthermore, a major limitation of several models is their dichotomic approach (slow and fast), yet there are experimentally measured Ca 2+ transients of at least four fiber types: I, IIA, IIX, and IIB [6]. Some models have not included important mechanisms dealing with Ca 2+ , such as the mitochondria or the NCX, despite their importance in shaping the Ca 2+ transients in different fiber types [12,13].
A recent model of skeletal muscle ECC, interestingly included the MITO and proteins such as the mitochondrial Ca 2+ uniporter (MCU) and the mitochondrial NCX (NCE) [14]. However, that model was performed based on Ca 2+ transients' measurements with Fura-2, while the most suitable dyes to study ECC in skeletal muscle seem to be the low affinity, fast Ca 2+ dyes, such as Mag-Fura-2 and Mag-Fluo-4 [1,15]. This explains why the reported amplitude of the signal (0.4 µM) was very low compared to previously published values being 10-20 µM in murine Flexor digitorum brevis (FDB) fibers with Mag-Fura-2 and Mag-Fluo-4 [15][16][17].
In this work, we used measurements carried out with the fast Ca 2+ dye Mag-Fluo-4 and a theoretical model to estimate the Ca 2+ movements produced during single and tetanic Ca 2+ transients in the four murine skeletal muscle fiber types. This permitted us to make a simulated comparison of the Ca 2+ movement across the four fiber types, including the variations in the mitochondrial [Ca 2+ ], the NCX, the SOCE, and their influence on the sarcoplasmic Ca 2+ regulation, and further overcoming the limitations imposed by slow Ca 2+ dyes.
Mathematical models that integrate information obtained on mammalian ECC under different experimental conditions have been used to simulate Ca 2+ transients in both slowtwitch and fast-twitch fibers [10,11]. However, the nature of the fibers used has not always been doubtlessly established. Furthermore, a major limitation of several models is their dichotomic approach (slow and fast), yet there are experimentally measured Ca 2+ transients of at least four fiber types: I, IIA, IIX, and IIB [6]. Some models have not included important mechanisms dealing with Ca 2+ , such as the mitochondria or the NCX, despite their importance in shaping the Ca 2+ transients in different fiber types [12,13].
A recent model of skeletal muscle ECC, interestingly included the MITO and proteins such as the mitochondrial Ca 2+ uniporter (MCU) and the mitochondrial NCX (NCE) [14]. However, that model was performed based on Ca 2+ transients' measurements with Fura-2, while the most suitable dyes to study ECC in skeletal muscle seem to be the low affinity, fast Ca 2+ dyes, such as Mag-Fura-2 and Mag-Fluo-4 [1,15]. This explains why the reported amplitude of the signal (0.4 µM) was very low compared to previously published values being 10-20 µM in murine Flexor digitorum brevis (FDB) fibers with Mag-Fura-2 and Mag-Fluo-4 [15][16][17].
In this work, we used measurements carried out with the fast Ca 2+ dye Mag-Fluo-4 and a theoretical model to estimate the Ca 2+ movements produced during single and tetanic Ca 2+ transients in the four murine skeletal muscle fiber types. This permitted us to make a simulated comparison of the Ca 2+ movement across the four fiber types, including the variations in the mitochondrial [Ca 2+ ], the NCX, the SOCE, and their influence on the sarcoplasmic Ca 2+ regulation, and further overcoming the limitations imposed by slow Ca 2+ dyes.

Results
Experimentally recorded, raw, single, and tetanic Mag-Fluo-4 Ca 2+ transients were calibrated in order to obtain the ∆[Ca 2+ ] in the sarcoplasm ( Figure 1). The [Ca 2+ ] peaks for the continuum of fiber types were: IIB and IIX: 16.58 µM, IIA: 12.77 µM, and I: 10.13 µM ( Figure 1A). For the tetanic Ca 2+ transients, subsequent peaks were also calculated (I: 11.24, 12.95, 14.81, and 16.47 µM; IIB: 12.13, 12.29, 12.24, and 11.96 µM) ( Figure 1B). ] during Ca 2+ transients elicited by a single and a train of AP (100 Hz) in mouse fibers at room temperature (21-23°C). All fibers were classified by myosin heavy chain direct determination in the original paper. Fluorescence recordings generated with Mag-Fluo-4 were used to obtain the ∆[Ca 2+ ] in the sarcoplasm produced by a single AP in fiber types I (green), IIA (blue), IIX (red), and IIB (black) (A) and a train of 5 AP in fibers I and IIB (B).
These calibrated signals fed all next estimations and simulations, as described in the Methods section and Table 1. First, we mathematically estimated the release rate of Ca 2+ (JRel) for both single and tetanic transients (Figure 2A,B; Table 2) and then simulated the Ca 2+ kinetics in the sarcoplasm ( Figure 2C,D), the SR ( Figure 2E,F), and MITO ( Figure  2G,H). Fibers type II peaked higher than fibers type I (~137% higher for IIA and 269-277% higher for IIX and IIB). The JRel estimated in tetanic Ca 2+ transients shows that the last peak's amplitude is reduced over 10 times for type I fibers and up to 7 times for IIB, IIX, and IIA ( Figure 2B; Table 2). The simulated sarcoplasmic ∆[Ca 2+ ] closely resembled the Figure 1. Experimental, calibrated measurements of [Ca 2+ ] during Ca 2+ transients elicited by a single and a train of AP (100 Hz) in mouse fibers at room temperature (21-23 • C). All fibers were classified by myosin heavy chain direct determination in the original paper. Fluorescence recordings generated with Mag-Fluo-4 were used to obtain the ∆[Ca 2+ ] in the sarcoplasm produced by a single AP in fiber types I (green), IIA (blue), IIX (red), and IIB (black) (A) and a train of 5 AP in fibers I and IIB (B).
These calibrated signals fed all next estimations and simulations, as described in the Methods section and Table 1. First, we mathematically estimated the release rate of Ca 2+ (J Rel ) for both single and tetanic transients (Figure 2A,B; Table 2) and then simulated the Ca 2+ kinetics in the sarcoplasm ( Figure 2C,D), the SR ( Figure 2E,F), and MITO ( Figure 2G,H). Fibers type II peaked higher than fibers type I (~137% higher for IIA and 269-277% higher for IIX and IIB). The J Rel estimated in tetanic Ca 2+ transients shows that the last peak's amplitude is reduced over 10 times for type I fibers and up to 7 times for IIB, IIX, and IIA ( Figure 2B; Table 2). The simulated sarcoplasmic ∆[Ca 2+ ] closely resembled the experimental recording described above. The [Ca 2+ ] SR rapidly decreased and slowly recovered as expected. Although qualitatively similar for all fiber types, quantitative differences arose mainly between the fibers type I and all fibers type II, in agreement with the fact that the Ca 2+ released by the fibers type I was the lowest. experimental recording described above. The [Ca 2+ ]SR rapidly decreased and slowly recovered as expected. Although qualitatively similar for all fiber types, quantitative differences arose mainly between the fibers type I and all fibers type II, in agreement with the fact that the Ca 2+ released by the fibers type I was the lowest. There was a higher value in ∆[Ca 2+ ]MITO for the fibers type I (up to four times) as compared to the three fibers type II. The total [Ca 2+ ] in MITO, considering both forms, bound to B and free, achieves a peak value of 0.26-1.12 µM after a single AP and 0.72-2.8 µM after 5 AP. Regarding the SR compartment, our calculations show that the available Ca 2+ is reduced up to ~63% (I: 82%, IIA: 63%, IIX: 64%, and IIB: 63%) for single transients, but to almost 53% (I: 76%, IIA: 63%, IIX: 56%, and IIB: 53%) after a train of five shocks ( Figure  2E,F).
We also calculated the variations in [Ca 2+ ] for the buffers present in the three simulated compartments (Figure 3), for both single (left column) and tetanic transients (right column). The differences in [CaPV] show the PV buffer influence in the first part of [Ca 2+ ] decay obtained in the IIX and IIB fibers, compared with I and IIA fibers. The second part of the decay, when the sarcoplasmic [Ca 2+ ] is near steady state, was less affected by the PV buffering. The Δ[CaDye] and Δ[CaATP] simulations were not shown as they resemble the shape of the sarcoplasmic [Ca 2+ ] already shown ( Figure 2C,D), but with different peaks (reported in Table 3). The ∆[CaDye] explains only 2.7% of the intracellular dye in the fiber type IIA and 2.2% in fiber types I, IIX, and IIB, thus ensuring dye unsaturation.   Figure S1, in order to show the return of [Ca 2+ ] to the basal levels.
There was a higher value in ∆[Ca 2+ ] MITO for the fibers type I (up to four times) as compared to the three fibers type II. The total [Ca 2+ ] in MITO, considering both forms, bound to B and free, achieves a peak value of 0.26-1.12 µM after a single AP and 0.72-2.8 µM after 5 AP. Regarding the SR compartment, our calculations show that the available Ca 2+ is reduced up to~63% (I: 82%, IIA: 63%, IIX: 64%, and IIB: 63%) for single transients, but to almost 53% (I: 76%, IIA: 63%, IIX: 56%, and IIB: 53%) after a train of five shocks ( Figure 2E,F).
We also calculated the variations in [Ca 2+ ] for the buffers present in the three simulated compartments (Figure 3), for both single (left column) and tetanic transients (right column). The differences in [CaPV] show the PV buffer influence in the first part of [Ca 2+ ] decay obtained in the IIX and IIB fibers, compared with I and IIA fibers. The second part of the decay, when the sarcoplasmic [Ca 2+ ] is near steady state, was less affected by the PV buffering. The ∆[CaDye] and ∆[CaATP] simulations were not shown as they resemble the shape of the sarcoplasmic [Ca 2+ ] already shown ( Figure 2C,D), but with different peaks (reported in Table 3). The ∆[CaDye] explains only 2.7% of the intracellular dye in the fiber type IIA and 2.2% in fiber types I, IIX, and IIB, thus ensuring dye unsaturation.  In general, the amount of Ca 2+ handled by the NCX was 41-121 times lower after 1 AP and 34-88 times lower after 5 AP, compared to the SERCA capacity. Moreover, the rate of transport by NCX was notably reduced when the sarcoplasmic [Ca 2+ ] achieves low values in IIX and IIB fibers given its low affinity. The SERCA pump maintains its influence in [Ca 2+ ] regulation near resting conditions. The [Ca 2+ ] returned to the sarcoplasm by the SOCE was negligible in all fibers compared to the total [Ca 2+ ] released. The total [Ca 2+ ] handled by each mechanism after the simulated time intervals, the free [Ca 2+ ] reached in each compartment, and the Ca 2+ bound to the chemical species are reported in Table 3. A longer time interval of the tetanic Ca 2+ transient simulation was included for some reactions ( Figure S1). The total amount of Ca 2+ obtained at rest considering all compartments of the model, including the extracellular space, remained constant during the simulated activation interval. We obtained that the variations in the total amount of Ca 2+ were lower than 10 −6 µM. This result evidence that truncation errors were negligible during the simulations.  In general, the amount of Ca 2+ handled by the NCX was 41-121 times lower after 1 AP and 34-88 times lower after 5 AP, compared to the SERCA capacity. Moreover, the rate of transport by NCX was notably reduced when the sarcoplasmic [Ca 2+ ] achieves low values in IIX and IIB fibers given its low affinity. The SERCA pump maintains its influence in [Ca 2+ ] regulation near resting conditions. The [Ca 2+ ] returned to the sarcoplasm by the SOCE was negligible in all fibers compared to the total [Ca 2+ ] released. The total [Ca 2+ ] handled by each mechanism after the simulated time intervals, the free [Ca 2+ ] reached in each compartment, and the Ca 2+ bound to the chemical species are reported in Table 3. A longer time interval of the tetanic Ca 2+ transient simulation was included for some reactions ( Figure S1). The total amount of Ca 2+ obtained at rest considering all compartments of the model, including the extracellular space, remained constant during the simulated activation interval. We obtained that the variations in the total amount of Ca 2+ were lower than 10 −6 µM. This result evidence that truncation errors were negligible during the simulations.

Discussion
The main findings of the present work were: (i) during a single twitch, the sarcoplasmic peak [Ca 2+ ] for fibers type IIB and IIX is between 15-25 µM, and for fibers type I and IIA reaches 6-12 µM, (ii) both the pattern of change and the peak concentrations of the Ca 2+bound species in the sarcomere, the sarcolemma, and inside the SR showed the order IIB ≥ IIX > IIA > I, (iii) the mitochondrial peak [Ca 2+ ] and the MITO buffers saturation showed the pattern I >> IIA >> IIX ≥ IIB.
Previous models of mammalian ECC were affected by either uncertainty in the classification of fiber types, unreliable kinetics of the raw Ca 2+ signals due to the use of slow Ca 2+ dyes, or the lack of information about the role of several intracellular compartments in Ca 2+ handling. To overcome these limitations, we based our model on the first calibration of Ca 2+ transients of the four main fiber types found in mammals, obtained using the fast Ca 2+ dye Mag-Fluo-4, and integrated new information gathered on MITO, SR, NCX, and SOCE, along basic knowledge on sarcoplasmic Ca 2+ movements and buffering.

A Model which Includes Four Fiber Types
Preceding models addressed ECC in one or two fiber types, mainly as most previous functional and biochemical information came from a dichotomic approach of muscle fibers: either slow vs. fast, or type I vs. type II [13,14,16,[34][35][36][37]. Additionally, simulations that have recently become spatially and mathematically more complex have remained biochemically oversimplified [38,39]. By taking Ca 2+ transients obtained in molecularly typed fibers covering the whole spectra from I to IIB as experimental source, our model adds novel information on the ECC-fiber type relationships.
The fact that molecular and biochemical differences (i.e., isoforms and their biochemical properties) underlie the differences in ECC among fiber types has been described in detail elsewhere and we refer the readers to those papers [1,6,40].
However, there is still a lack of molecular and biochemical information, particularly for fibers type IIA and IIX. Our adjustments and results for these two types of fibers seem reliable as most values obtained laid between those of fibers type I and IIB. This agrees with the fact that fibers IIA and IIX showed kinetics of the Ca 2+ transients, Ca 2+ sensitivity, and other dynamic properties which are mostly intermediate between I and IIB [6,[40][41][42]. Together, this confirms that most values of molecular, biochemical, and physiological parameters follow a continuum from I, to IIA, to IIX, to IIB.

[Ca 2+ ] Kinetics with a Fast Ca 2+ Dye
The Ca 2+ transients modeled in the present study were obtained with the fast Ca 2+ dye Mag-Fluo-4. Two issues with quantitative impact on the results deserve attention: the biochemical properties of the dye and its loading conditions. This dye has a 2:1 Mag-Fluo-4-Ca 2+ binding stoichiometry, with an in situ K d of 1.652 × 10 5 µM 2 [17]. This explains why it can reliably track the Ca 2+ transients even in the fastest types, i.e., IIB, demonstrate subtle differences among all four fiber types and resolve every peak in a tetanic transient, being a trustable source for the model. On the other hand, as the fibers typically have~200 µM of Mag-Fluo-4 [17], we found that less than~3% of the dye is bound to Ca 2+ .
The very low affinity of the dye and the lack of saturation confer Mag-Fluo-4 the ability to determine a trustable peak sarcoplasmic [Ca 2+ ]. Our results with Mag-Fluo-4 (calibration performed in the present study and [17]) join those with Mag-Fura-2 [16,35] which show that the peak [Ca 2+ ] in fast fibers IIX or IIB is typically between 15 and 25 µM in mammalian skeletal muscle (16-22 • C). The numbers for type I and IIA fibers also agree with a study reporting peak Ca 2+ for soleus slow fibers [13]. That work, however, may have really used type I and IIA fibers, as the periphery of the soleus muscles of mouse used in the experiments have both types of fibers evenly distributed, as we have verified using specific antibodies and confocal microscopy (not shown). In conclusion, that paper and our results agree in that peak [Ca 2+ ] for fibers types I and IIA is between 6 and 12 µM [13]. The above numbers are one order of magnitude higher than those reported with slow dyes [14,33,43,44]. This fact may, for instance, reduce the estimated maximum rate of Ca 2+ flux during SOCE or excitation-coupled Ca 2+ entry (ECCE) activations [33,45], as this value depends on the driving force for Ca 2+ [46]

Comprehensive Integration of Mechanisms Involved in Ca 2+ Handling: Sarcoplasm, SR, MITO, NCX and SOCE
The ECC is becoming more complex [1] and experimentally addressing some questions is challenging. For instance, important papers have investigated Ca 2+ kinetics in MITO and SR, as well as SOCE function in fast fibers [33,47,48], but differences among all fiber types are infrequently [12] or never studied. The mathematical model presented here helps address these limitations.
The peak flux differences found in the release rate are consistent with the two-to ten-fold higher content of the ryanodine receptor (RyR) in the fastest fibers, compared to the other fibers [6,[49][50][51]. The comparable values of half-width and rise time in all fiber types can be explained as all fiber types share the RyR1 isoform [52].
Upon release, Ca 2+ has multiple fates: (i) binds to classical sarcoplasmic buffers such as PV, Tn, TNS, ATP, and the dye; (ii) enters into MITO and SR, or (iii) is recycled through the NCX and SOCE mechanisms [1].
In fast fibers (16-22 • C), [CaTn] has ranged from 80 µM [27] to around 240 µM [16,39], being~199 µM for IIA, IIX, and IIB in the present study, within 7 ms from the beginning of the release of Ca 2+ . Slow fibers gave values of 85 µM [13] and~110 µM (fibers type I in the present study) within 15 ms. Different reaction rate constants arising from temperature corrections, our larger SERCA capacity, as well as the inclusion of MITO and SOCE, explain why the peak [CaTn] is somewhat lower in the fast fibers in our model than in the model of Baylor and Hollingworth (2007) [16]. The TNS peaked 31 µM during a train of APs in slow fiber types after 200 ms [13], and in the present work, such as with Tn, have peaked below, being 22 µM. Specifically in type IIX and IIB fibers, which have a higher quantity of PV, their influence is reduced, as the binding and unbinding rates are similar to those of PV.
The peak concentration of [CaPV] in fast fibers has ranged from 90 µM [27] to 120 µM [16], close to the 116-142 µM for IIX and IIB found in the present study, after 25 ms of Ca 2+ release. Slow fibers are devoid of PV [13], or very low values have been measured in fibers type I and IIA in murines [53,54], justifying the values found (~0.5-5 µM) here. Then, the continuum of [CaPV] in I, IIA, IIX, and IIB fibers followed the order IIB > IIX >> IIA > I.
The peak [CaATP] we found for fibers IIX and IIB coincides with the values reported previously for fast fibers, i.e., 60 µM [27] and~70 µM [16]. In all cases, the peak was within 6 ms from the release. In type I fibers, we report a value of 23 [14]. The peak [Ca 2+ ] was~0.25 µM in our results for fibers type IIX and IIB, about 25 ms after Ca 2+ release. This concentration is similar to that obtained experimentally in fast mouse fibers [47]. Somewhat later, [Ca 2+ ] peaked inside MITO in type I fibers, reaching a value of 0.9 µM. This delay likely reflects Ca 2+ diffusion into MITO, as measured in mouse using genetically encoded sensors [55]. Notably, the peak [Ca 2+ ] was far higher in I than in II, likely as the MCU of slow fibers has more capacity than that of the fast fibers [30], with no difference in the properties of the NCE, allowing more Ca 2+ to accumulate inside MITO of fibers type I. Additionally, in fibers type I, the lower amount of PV compared to II allows more Ca 2+ free to diffuse into MITO. For fibers IIA, their intermediate amount of Ca 2+ released, along with their intermediate amount of PV and SERCA content may explain their midsize value of Ca 2+ inside MITO. According to their peak amplitude, Ca 2+ transients in MITO can be ordered as I >> IIA >> IIX ≥ IIB.
Previous works have estimated the total concentration of MITO buffers [B]. From the simulated analysis of the total [B] in cardiac muscle fibers performed in [56], a value of 2 µM was found to be optimal. In the work of Marcucci et al. (2018) [14] for skeletal muscle fibers a larger buffer concentration of 20 µM was also tested. In the present study, a value of 20 µM was used. A sizeable total [B] inside MITO is expected in order to deal with the larger and faster amounts of Ca 2+ released in skeletal compared to cardiac muscle, and to explain the slowing of decay of the Ca 2+ transients when the MITO uptake is blunted [12,47]. However, well calibrated, high resolution, experimental measurements of MITO Ca 2+ buffering and kinetics are still required to reach an accurate estimate of the total [B] and [Ca 2+ ] inside MITO in skeletal muscle fibers. Ca 2+ pumped by the SERCA was estimated to be 1.5-3.5 µM in fast fibers [16], but amounted up to~50 µM for IIX fibers in our model, after 25 ms. In slow fibers, a value somewhat below 1 µM was reported [13] at the same time (50 ms), but we found~15 µM for type I. The differences between models can be explained due to temperature, as well as half-width differences in the recordings that fed the simulation. According to the intermediate kinetics of the decay phase of the Ca 2+ transients in fibers IIA, as well as our previous discussion, midsize values for IIA were expected. The difference of Ca 2+ pumped by SERCA 50 ms after Ca 2+ release was IIB = IIX >> IIA >> I. The total capacity of Ca 2+ extrusion by the NCX in fibers type IIB, IIX, and IIA was only 1.1-1.3 times higher than in I, suggesting that the larger amount of NCX1 in fibers type I is balanced by the larger capacity of the NCX3 present in fibers type II. The Ca 2+ extruded by the t-system seems to be immediately recycled and has been called a counter-flux [5]. Our results confirmed that SOCE is quantitatively small in skeletal muscle [33,46], and it is difficult to speculate on the importance of the minor differences among fibers. More robust experimental data should be gathered before a conclusion about this issue can be stated.

Final Remarks
Although the above analyses give averages of peak sarcoplasmic/compartments [Ca 2+ ], spatially refined models for fast fibers have shown up to a 20-fold gradient in the sarcoplasmic [Ca 2+ ], depending on the distance of a subcellular region from the Ca 2+ release units [13,16,39]. This phenomenon is also expected to apply to all fiber types, but the magnitude of those gradients inside the fibers was not explored in the present study.
Single-compartment simulations are expected to have errors associated with an inability to estimate local gradients in [Ca 2+ ] and in the [Ca 2+ ] bound to the binding sites. In the present study, the SOCE flux is associated with changes in [Ca 2+ ] in the SR compartment. However, Ca 2+ gradients in the SR have been measured during ECC and the Ca 2+ levels decrease more rapidly in the terminal cisternae than in other regions of the SR [57]. Therefore, a model that considers the gradients in the SR would allow a more accurate estimate of the amount of Ca 2+ entered by the SOCE.
Additionally, as our Ca 2+ -bound chemical species had higher concentrations than those recently estimated, the thermal changes associated with ECC in mammalian muscle should be higher than proposed [27]. Our model may be a source to build a more complete model on thermal changes in all fiber types during single and tetanic stimulation.

Experimental Single and Tetanic Ca 2+ Transients
Typical fluorescence experimental recordings (F) of single Ca 2+ transients of fiber types I, IIA, IIX, and IIB were taken from [6]. Tetanic Ca 2+ transients (100 Hz) of types I and IIB fibers were taken from [12]. The tetanic Ca 2+ transient of IIA and IIX fibers were simulated assuming that they share morphology with I and IIB fibers, respectively, as published [6,12]. In all cases the signals were obtained with Mag-Fluo-4 in electrically stimulated isolated fibers from mouse muscles. Fibers were classified based on myosin heavy chain determination as detailed described in the above references. The conversion of F of single and tetanic experimental recordings to [Ca 2+ ] was performed according to the calibration method and the value parameters presented in [17]. Briefly, [Ca 2+ ] was calculated using the expression: where F max and F min are the maximum and minimum fluorescences (150.9 A.U and 0.14 A.U) respectively, K d is the in situ dissociation constant of Mag-Fluo-4 (1.652 × 10 5 µM 2 ), and [D] T is the total dye concentration (229.1 µM), at 20 • C. As the experimental Ca 2+ transients that feed the model were acquired between 21-23 • C, no temperature corrections were performed in the calibration parameters.

Model Description
The rate of change of free [Ca 2+ ] was described in three compartments (SR, sarcoplasm, and MITO) with a system of differential equations. The rate of change of [Ca 2+ ] with respect to time in the SR was determined by where J Rel is the release rate flux, J SERCA is the SERCA flux, and F the reaction rate of Ca 2+ in the SR with calsequestrin (CSQ). The rate of change of sarcoplasmic [Ca 2+ ] was expressed as: where J MCU was the inflow through the MCU, J NCE outflow through the NCE, and F([Ca 2+ ] MITO , B) the reaction rate with B. The B was 20 µM of binding sites, which react with Ca 2+ [14]. The parameters that describe the J Rel , and the values of maximum capacity of SERCA, NCX, NCE, and SOCE (V SERCA , V NCX , V NCE , and V SOCE ), were fitted to the measured Ca 2+ transients. All other values of the model were taken from the literature, and the specific references are given in Table 1 and Table S1.

Release Rate of Ca 2+
The J Rel of Ca 2+ from the SR in fast and slow fibers were estimated from measurements of sarcoplasmic Ca 2+ transients obtained with Mag-Fluo-4 ( Figure 1A), as previously done with Mag-Fura-2 [35]. This approximation assumes that Ca 2+ and other chemical species that react with Ca 2+ are uniformly distributed in the sarcoplasm, which is thus described as a single compartment. The time derivative of the total Ca 2+ concentration in the sarcoplasm, d[Ca] T /dt, is given by the sum of the inflows and outflows through the membrane. Therefore, from the sum of d[Ca] T /dt and outflows, we obtained the release rate of Ca 2+ . A gaussian model was used to simulate the J Rel elicited by a single AP: where N is the number of peaks used to fit the J Rel produced by a single AP, R is the peak amplitude, T i is the location in the time axis, and τ relates to the width of each peak. R, T, and τ were adjusted to fit the peak amplitude, half-width, rise time, and decay time of the release rate of Ca 2+ elicited by a single AP (Table 2). f Rel is the multiplication factor used in the second and subsequent APs to fit the measured peaks amplitudes. M is the number of APs and T 2 of 10 ms is the time between stimulations.

Reaction of Ca 2+ with Sarcomeric Buffers
The local reactions of Ca 2+ with S of the chemical species in the compartments were described by the law of mass action,  [14]. A temperature of 22 • C and a Q 10 of 2 were considered to adjust the reaction rates (Table S1). The resulting ordinary differential equations system was solved using the ode15s solver in MATLAB 2021b (MathWorks, Natick, MA, USA). The occupancy fraction of the binding sited at equilibrium were obtained from the simulated data.

Muscle Proteins Concentration
The Tn molecules concentration were 120 µM in all fiber types such as in Hollingworth et al. (2012) [13]. As each fast fiber Tn molecule has two Ca 2+ binding sites and slow fiber molecules one, we thus assumed a concentration of 240 µM with positive cooperativity for IIA, IIX, and IIB fibers, and 120 µM for type I fibers. Additionally, as each Tn molecule has two non-specific sites (TNS) we considered concentrations of 240 µM for the binding sites of fast and slow fibers. As the Tn molecules are located in the myofibrillar space (MS), the average [Tn] and [TNS] in the sarcoplasm was rescaled by the occupation of the MS volume in the sarcoplasm (V MS and V sarc in Table 1).
The [PV] in the fastest fibers can achieve large concentrations, 1.5-1.9 mM [16,27]. We also considered that fibers type I have 300 times less [PV] [7], leading to 6 µM. Fibers IIA have a [PV] closer to fibers type I than type IIB, so we assumed 10 times more PV content in IIA than in I type [58], giving 60 µM, which is far lower than IIX and IIB, as expected. The ATP, Dye, CSQ, and B concentrations were listed in Table 1.

Reuptake Rate of Ca 2+ by SERCA
The reuptake of Ca 2+ by SERCA to the SR is described by the Michaelis-Menten kinetics and expressed as: where V SERCA is the maximum flux rate, [Ca 2+ ] is the Ca 2+ concentration in the sarcoplasm, K SERCA is the dissociation constant and h SERCA the Hill coefficient. The SERCA isoform predominantly expressed in fast-twitch muscle fibers is the 1a, whereas in slow-twitch is the 2a [59]. Both have similar K SERCA and h SERCA [29], however different content of SR pump molecules [8], and thus V SERCA . The V SERCA was then adjusted to match the decay phase of the measured and simulated Ca 2+ transients in all fiber types.

The Mitochondrial Ca 2+ Uniporter (MCU) Inflow and the Sodium Ca 2+ Exchanger (NCE) Outflow
The J MCU was represented by a saturable first-order transporter, independent of the internal [Ca 2+ ] MITO : (8) where V MCU is the maximum flux rate, [Ca 2+ ] is the Ca 2+ concentration in the sarcoplasm, and K MCU is the [Ca 2+ ] where the transport rate is half-maximum. V MCU values were taken as 0.012 µM ms −1 in fast and 0.049 µM ms −1 in slow fibers at 16 • C [13,16], and adjusted considering a Q 10 of 2.
The J NCE was modelled with a stoichiometry of 3:1 and described with the following expression [31]: where ∆Ψ m,mito is the mitochondrial membrane potential, V NCE is the NCE activity, and K NCE,Ca 2+ and K NCE,Na + are the dissociation constants for the Ca 2+ and Na 2+ binding to the NCE, respectively. F, R, and T denote the Faraday constant, the ideal gas constant and temperature, respectively. V NCE was modulated in the fastest fibers to reproduce the speed of the [Ca 2+ ] mito decay phase, completed in a period of about 100 ms and measured during a single twitch, as in Rudolf et al. (2004) [55]. 4.2.6. Ca 2+ Flux through the NCX To describe the sarcolemmal NCX, the same expression of the NCE, with stoichiometry of 3:1 was used, although with different simulated values. Reported Ca 2+ -transport rate values for SERCA range from 10 (in membrane vesicles) to 13-14 nmol mg −1 min −1 (NCX3 and NCX1 isoforms, respectively) [60,61]. We assumed a maximum transfer rate (V NCX ) between 30 and 40% higher than that of SERCA. As the NCX is predominantly located in the T-tubule membrane [62], whereas the SERCA in the SR membrane, we thus rescaled V NCX by a ratio (T-tubule membrane/SR membrane) of 0.066 and 0.111 in fast and slow fibers [20], i.e., between 214.5 and 120.1 µM s −1 for all fibers type II. In fibers type I we used 93.2 µM ms −1 . For K NCX,Ca 2+ and K NCX,Na + , we, respectively took values of 140 µM and 14 mM in fast fibers, and 130 µM and 11 mM in slow fibers [32].

The Effect of Store-operated Ca 2+ Entry (SOCE)
We described the flux, J SOCE , as a fraction, P SOCE , of a given maximum value, and V SOCE as J SOCE = V SOCE P SOCE [63]. As said, P SOCE change as a function of the [Ca 2+ ] in the SR, which can be described by the function: where K SOCE and h SOCE are the Ca 2+ dissociation constant of STIM1 (0.35 mM) and the Hill coefficient (4.7), respectively [33]. We assumed that the V SOCE in the fastest IIB fibers achieved up to 35 µM s −1 [33]. Lower values were used for IIX, IIA, and I fibers.

Conclusions
Our mathematical, comprehensive model allows us to gain insight into the kinetics of the Ca 2+ transients obtained with the fast Ca 2+ dye Mag-Fluo-4, for the continuum of skeletal muscle fiber types. Sarcoplasmic peak [Ca 2+ ] is one order of magnitude higher than reported with slow dyes. The magnitudes of change of the Ca 2+ -bound forms of the Ca 2+ buffers studied follow the order IIB ≥ IIX > IIA > I, except for mitochondrial peak [Ca 2+ ] which showed the pattern I >> IIA >> IIX ≥ IIB. The kinetics for fibers IIA and IIX proved to be intermediate between I and IIB fibers, supporting dynamic data. The results may help better quantitate SOCE fluxes and thermal changes in mammalian fiber types in the future and support the use of fast Ca 2+ dyes for most experimental approaches in skeletal muscle.