Incorporation of sarcolemmal calcium transporters into the Shorten et al. (2007) model of skeletal muscle: equations, coding and stability

RETROSPECTIVE Abstract We describe a major development of the Shorten et al. (Shorten et al., 2007) model of skeletal muscle electrophysiology, biochemistry and mechanics. The model was developed by incorporating equations for sarcolemmal transport of calcium ions, including L-type calcium channel, sodium-calcium exchange, calcium pump and background calcium channel. The extended model also includes an addition to the equations for extracellular potassium ion movements to enable the exchange of potassium ions between bulk (plasma) concentration and the interstitial and tubular compartments to be modelled. In further research in an accompanying paper (Tasaki et al, 2019), we succeeded in reproducing muscle cramp, as well as its prevention and reversal, by investigating muscle contraction and cramp using this extended model in comparison with the original model.


Introduction
The aim of this paper is to extend the Shorten et al. skeletal muscle model (Shorten et al., 2007) to include surface membrane calcium transporters. This is necessary to enable the model to be used in applications involving interactions between the surface transporters and either drugs or ions. The extended model was developed specially for a project involving a medication that acts on L-type calcium channels described in an accompanying article (Tasaki et al. 2019).  build-up (membrane) and no Pi feedback on the cross-bridge dynamics (membrane+precipitation). Also shown are the calcium and Pi dynamics in a fibre exhibiting total fatigue response (membrane+precipitation+cross-bridge). Once the Pi exceeds a critical threshold it precipitates with SR calcium to reduce calcium release from the SR and cross-bridge cycling kinetics. B -the measured EDL response to sustained 25 Hz stimulation.

Original model
The skeletal muscle model developed by Shorten et al. (Shorten et al., 2007) is extensive. It includes models of the action potential, electromechanical coupling, the cross-bridge reactions, and the immediate biochemical processes of consumption of ATP and concentrations of inorganic phosphate. The model was also parameterised to fast and slow twitch fibre types, which have a number of biochemical and biophysical differences. These differences are related to the different force/fatigue responses in fast and slow twitch fibres that are not well understood. A better understanding of the cascade of biochemical and physical events underlying skeletal muscle fatigue will provide a basis for targeted intervention by nutritional or other means to improve muscle function. Loss of muscle function is a major issue in elite athletes, cachectic subjects and those undergoing aging, where muscle function decreases after the age of 40.
The model is available in full working order from the Physiome Model Repository (PMR) at https://models.physiomeproject.org/exposure/487ba7abea1b0ed7ee2ee1bef3f13aee Since the model includes the immediate biochemical energetic processes, including variations in ATP and phosphate, it can reproduce the sequence of events underlying muscle fatigue, as 2/11 shown in Figure 1. Phosphate is generated during cross-bridge cycling.

Explanation why sarcolemmal calcium transporters need to be incorporated
Electromechanical coupling in smooth and cardiac muscle depends on calcium influx. In both types of muscle, release of calcium ions from the sarcoplasmic reticulum is triggered by calciuminduced calcium release. By this process, a small calcium influx is amplified to produce sufficient calcium to initiate cross-bridge reactions and so produce sliding filament contraction. By contrast, electromechanical coupling in skeletal muscle is purely electrical. Tight coupling between T-tubule membranes through which the action potential propagates makes it possible for SR calcium release to occur even if no calcium influx occurs through the sarcolemmal membrane. Skeletal muscle contractions can therefore continue to occur even in the absence of extracellular calcium. This was demonstrated by Edman and Grieve (Edman and Grieve, 1964). The immediate effect of exposure of a frog skeletal muscle to nearly zero external calcium was that the muscle continued to contract. Cairns et al. (Cairns et al., 1998) found a similar result in mammalian skeletal muscle. Nominally, calcium-free solutions "had little effect on tetanic force in non-fatigued muscle." Figure 2. Effect of lowering calcium concentration in the solution bathing a frog sartorius muscle on mean resting potentials (upper curves; total range of potentials also shown), tetanic tensions (middle curves) and the twitch:tetanus ratio (lower curves). Frog sartorius muscle, temperature 0-1 o C. Square pulses of 1 ms duration were applied through two wire electrodes so as to stimulate the whole muscle. Calcium-free solution contained 0.1 mM EDTA. (Reprinted with permission from Edman and Grieve (Edman and Grieve, 1964)).
However, during prolonged exposure to very low calcium, the twitch contractions fall. In Figure 2, the muscles were exposed to either 0.1 mM, 0.01 mM calcium or calcium-free Ringer's solution.
Both tetanic tension and the twitch:tetanus ratio slowly declined. This result implies that over long periods of time skeletal muscle does depend on sarcolemmal calcium influx. See Cairns (Cairns et al., 2003;Cairns and Lindinger, 2008) for more discussion on the role of extracellular calcium.
It is therefore important to incorporate surface membrane calcium fluxes in the Shorten et al. model.

Methods of incorporating sarcolemmal calcium transporters L-type calcium current
The CaV 1.1 and CaV 1.2 L-type calcium channels are found in adult skeletal muscle (Jeftinija et al., 2007;Flucher and Tuluc, 2017) . CaV 1.1 is restricted to the T-tubules and CaV 1.2 is expressed in the sarcolemma of type I and IIa myofibers, but not IIb fibres. The data and equations for the L-type calcium current (ICaL) were obtained from the work of Wang et al. (Wang et al., 2008) who studied the effect of a sodium channel blocker, Riluzole, on the action potential in cultured human skeletal muscle cells. The equations were taken from the Physiome Model Repository: https://models.physiomeproject.org/exposure/3a4c502f78a8c3a5845844b501e153e9 The equations are: (1) The equation for the inactivation (f i nf ) process was added for use in a future project.

Sodium-calcium exchange
The data and equations for sodium-calcium exchange (NCX) were obtained from Noble et al. (Noble et al., 1991), which were originally developed by DiFrancesco Noble (Di Francesco and Noble, 1985) and which reproduce the data obtained by Kimura, Miyamae and Noma (Kimura et al., 1987) for cardiac muscle. NCX is based on 3:1 stoichiometry (Michel et al., 2014). NCX is also more abundant in slow twitch fibres, which is consistent with a role for NCX in fatigue resistance (Michel et al., 2014).

Calcium pump
The data and equations for the calcium pump current, IpCa, were obtained from Lindblad et al. (Lindblad et al., 1996) in cardiac muscle. The equation is:

Background calcium leak
The data and equations for the background calcium current leak, IbCa, were obtained from Noble et al. ((Noble et al., 1991)).
The equation is: g bC a was usually set to zero. The calcium pump and sodium-calcium exchange were sufficient to balance calcium movements in steady state conditions.
To represent the bulk potassium concentration outside the muscle, [K]b (which is equivalent to a plasma potassium concentration), we added a term to the differential equation for interstitial potassium, K e , which then becomes: [K] b is assumed to be fixed.
The original CellML file along with all the codes can be found in the following link in the PMR: https://models.physiomeproject.org/workspace/5c6 It is worth to mention that in order to be able to run python scripts the corresponding sedml and cellml files need to be downloaded in the same folder. Also In order to be able to run sedml file the corresponding cellml file needs to be downloaded in the same folder.

5/11 3 Model results
All simulations were run using OpenCOR (version 0.6) (Garny and Hunter, 2015) (http://opencor.ws) Comparison with the original model over short pulse trains The curves also compare well qualitatively with Figures 6 and 11 in the original model paper. There is, however, an important quantitative difference arising from the fact that we used initial conditions obtained from running the extended model for very long periods (hours) to achieve steady state conditions. With the original model parameters, this produces a lower initial state of intracellular calcium and, therefore, much lower contraction than in the Shorten et al. paper.
We also ran these computations using the initial conditions, i.e. not steady state conditions, obtained from the CellML file for the original model on PMR. These results were very similar qualitatively to those shown in Figure 3, except that the contraction was much larger. There was 6/11 nevertheless still a difference from the published Shorten et al. paper. The contraction was found to be about 50% of that in the original paper.

Test of calcium balance
Although the equations for influx of calcium through ICaL were obtained on skeletal muscle, the equations for efflux of calcium via sodium-calcium exchange and the calcium pump were obtained from heart muscle models. NCX1 is expressed in high levels in cardiac muscle. NCX1 is also found in skeletal muscle (Cifuentes et al., 2000), where it is found in sarcolemma and T-tubule membranes (Sacchetto et al., 1996). The major function of the skeletal muscle NCX appears to be to transport calcium out of the cell after contraction.
Their role in the model is to balance the influx and efflux so that over many cycles of normal physiological muscle activity there would be no net change in calcium content of the muscle fibres. We adjusted their amplitudes until balance was achieved, as shown in Figure 4. The simulation file Fig04.sedml contains the computational setting for running the model.  Figure 5 shows that, when the L-type calcium channel is blocked, the model is no longer in a steady state. SR calcium, free intracellular calcium and contraction all decline slowly, as observed by Edman and Grieve (Edman and Grieve, 1964) in their experiments in skeletal muscle in low or zero external calcium. The model achieves a new steady state after 1-2 hours. This very slow change and its experimental validation is dealt with in a further paper (Tasaki et al., 2019). The simulation file Fig05.sedml contains the computational setting for running the model. In addition to incorporating equations for surface membrane calcium fluxes, we added a term to the extracel- lular potassium equations to represent the movement of potassium ions between extracellular spaces (interstitial and tubular) within the muscle structure and the spaces outside the muscle, which would include blood plasma potassium. This feature of the extended model is used in an accompanying paper (Tasaki et al., 2019) to model the effects of potassium wash-out during increased blood flow. Figure 6 shows that the equations behave as expected. The bulk (plasma) potassium, [K]b, was switched suddenly between 5 and 2 mM twice during the calculation. The time constant, TauK3, was arbitrarily set to 4 seconds to represent fast perfusion of the interstitial space. Interstitial potassium correctly follows an exponential decline when [K]b is switched to 2 mM and recovers equally rapidly when [K]b is returned to 5mM. The use of this process in representing perfusion of the interstitial space is shown in the accompanying paper (Tasaki et al., 2019). The simulation file Fig06.sedml contains the computational setting for running the model.

Deficiency of the extended model
As already noted, when the results shown in Figures 4 and 5 are examined carefully, the magnitude of the contractions is clearly very much smaller than in the original Shorten et al. model. The computed steady state cross-bridge reactions are much less than 1 µM, whereas in the original model the reactions are at least an order of magnitude larger. The reason for this discrepancy is that the model was not parameterised for very long runs, of hours rather than minutes. All our computations are performed when the extended model has achieved steady state after about 2 hours since we need a steady state condition to study the dynamics of change from the steady state when the model is perturbed, e.g. by partial or complete block of the L-type calcium channel, as done in the accompanying paper (Tasaki et al., 2019).
We therefore investigated whether a simple parameter change could restore the magnitude of steady state contraction after very long time periods. Since the low value of calcium and contraction after 2 hours represents a loss of sarcoplasmic reticulum calcium, we investigated the effect of increasing the uptake rate for calcium entry into the SR, nuSR. Figure 7 shows that this simple parameter change can indeed increase the steady state contraction by more than an 8/11 However, the results also showed that further increase of nuSR (e.g to 45) prevented the membrane potential from fully repolarizing. The reason is that the increased calcium release during each contraction activates the sodium-calcium exchange, so producing a larger after-potential.
There is therefore a deficiency that requires further re-parameterisation of the extended model in future applications.

Discussion
Sarcolemmal calcium influx is not necessary for individual twitch responses in skeletal muscle since the membrane trigger for electromechanical coupling is entirely electrical, independent of calcium current. This property is what enabled Shorten et al. (Shorten et al., 2007) to reproduce fatigue in their model with only intracellular recycling of calcium.
However, by adding equations for the sarcolemmal calcium fluxes we have shown that over periods of time greater than a few seconds, the sarcolemmal calcium fluxes would be expected to have a significant effect. Specifically, the small fluxes gradually alter SR calcium and intracellular free calcium. These changes in calcium concentrations then produce changes in contraction. The slow time course of those changes corresponds to the dynamics of change in contraction observed experimentally when calcium influx is reduced, as shown in the accompanying paper (Tasaki et al., 2019). The influences of changes in sodium, potassium and calcium ion concentrations, on the contraction of skeletal muscle are very complex (Cairns et al., 1998(Cairns et al., , 2003Cairns and Lindinger, Figure 7. Computed relation between the SR uptake rate, nuSR, and the steady state contraction at a frequency of stimulation of 0. Further work also needs to be done to explore parameter sets that would increase the steady state contraction to be comparable to the contraction level calculated in the original model. A possible way forward to be developed in the future would be to rebalance the proportion of calcium efflux attributable to sodium-calcium exchange compared to the sarcolemmal calcium pump. The sodium-calcium exchange produces a depolarizing current, whereas the calcium pump does not do so. We have not systematically explored this question in this paper since it is likely that other parameter adjustments may also be required. The parameter space to be explored may be quite large. There are details/documentation on how the source code was compiled There are details on how to run the code in the provided documentation The initial conditions are provided for each of the simulations Details for creating reported graphical results from the simulation results Source code: a declarative language is used (e.g. SBML, CellML, NeuroML) The algorithms used are defined or cited in previous articles The algorithm parameters are defined Post-processing of the results are described in sufficient detail

Executable model provided:
The model is executable without source (e.g. desktop application, compiled code, online service) There are sufficient details to repeat the required simulation experiments

The model is described mathematically in the article(s):
Equations representing the biological system There are tables or lists of parameter values There are tables or lists of initial conditions

Machine-readable tables of initial conditions
The simulation experiments using the model are described mathematically in the article: Integration algorithms used are defined Stochastic algorithms used are defined Random number generator algorithms used are defined Parameter fitting algorithms are defined The paper indicates how the algorithms yield the desired output