OPEN ACCESS Reproducible Model Model of skeletal muscle cramp and its reversal

In an other paper (Noble et al., 2019), we developed the Shorten Shorten et al. (2007) model of skeletal muscle by incorporating equations such as surface calcium ﬂuxes. In further research in this paper, we succeeded in reproducing muscle cramp, as well as its prevention and reversal, by investigating muscle contraction and cramp, in which calcium regulatory networks are involved, using the extended model in comparison with the original model. Incorporation of data from a Traditional Medicine from root extracts of peony and licorice and one of its pure chemicals was modelled. Sensitivity analysis of the extended model shows the robustness of the calcium regulatory networks. Muscle cramp, in the extended model, requires calcium inﬂux via the L-type calcium channel and it will not occur without calcium inﬂux. Reduced calcium inﬂux can delay or prevent cramp. Increased interstitial potassium is implicated in developing and maintaining cramp. Mechanism of reversal of cramp requires wash-out of extracellular potassium via increased blood ﬂow, followed by calcium eﬄux via sodium-calcium exchange. This paper shows the ﬁrst successful quantitative electrophysiological and mechanical model of cramp and of its reversal. 2+ eﬄux, Na + - Ca 2+ exchanger.


Purpose of development of the original model
Both fatigue and cramp are important pathologies of skeletal muscle. The Shorten et alShorten et al. (2007) model reproduces muscle fatigue, but not muscle cramp. Yet, the latter is incapacitating and often painful. Without such a model, it is impossible to make or test predictions concerning the prevention or reversal of muscle cramp.
The prevalence of muscle cramp in patients with muscle cramp over the age of 60 is around 46% (Maisonneuve et al., 2016), with a high impact on quality of life: 21% of sufferers "described their symptoms as very distressing". It is therefore a widespread health problem in aging populations. This data is for rest cramps (Naylor and Young, 1994). Exercise associated muscle cramps are also common with acute pain, and distress affecting performance that can last days (Miller et al., 2010). The authors note that "their cause remains unknown".
It is difficult to account for the development of muscle cramp without supposing that calcium influx may be important. Cramp may be defined as a maintained contraction that cannot be immediately reversed by ceasing to electrically excite the muscle. Overloading the muscle fibres with calcium might satisfy this criterion since it would maintain contraction for as long as the free intracellular calcium is above the threshold for activating the cross-bridge reactions of the contractile proteins. Figure 1. Actions of SKT in isolated diaphragm of guinea pigs in twitch amplitude at a frequency stimulation of 0.1 Hz. A: top trace:1 mg/ml SKT shows no action; bottom trace: 10 mg/ml SKT shows slow action. B: dose-response relation of actions of SKT in isolated diaphragm of guinea pigs in twitch amplitude at a frequency stimulation of 0.1 Hz (from ) It is therefore important to incorporate surface membrane calcium fluxes in the Shorten model (Shorten et al., 2007). As we have shown in the other paper (Noble et al., 2019), our extensions of the model succeed in incorporating these fluxes, as well as those of extracellular potassium ions in exchange with bulk (plasma) potassium, [K] b .

Supporting experimental data showing why sarcolemmal calcium transporters need to be incorporated
For interaction with this modelling research, we use some of the experimental data performed in the Oxford laboratories on a 1:1 combination of root extracts from peony and licorice and one of its pure chemicals, glycyrrhetinic acid (GA). Its formula first appears in the Shang-han-lun (CE c 200) for muscle disorder and formulae of peony and licorice appears in the Ishinho (984). In Japan, it is called Shakuyakukanzoto (SKT), which has been used for relieving muscle cramp ( (Hyodo et al., 2006)). From the following data, we hypothesize that reduction of Ca 2+ influx is one of the important actions in reducing muscle contraction. First, in isolated skeletal muscle (diaphragm of guinea-pig), the twitch contraction in response to electrical stimulation, had a very slow reduction at high doses of SKT (either through nerves or directly applied to the muscle) (Figure 1 A). We calculate that the time constant of decline in contraction by 10 mg/ml SKT is around 21 minutes ( ). The dose-response relation of actions of SKT in isolated diaphragm of guinea pigs in twitch amplitude at a frequency stimulation of 0.1 Hz shows that the relation is steep (Figure 1 B)( ) Second, in experiments on the pore-forming subunit of the skeletal muscle L-type calcium channel expressed in Human Embryonic Kidney (HEK) cells, fluorescence imaging experiments show that relative calcium influx is reduced by SKT (Figure 2 A) and by one of its pure chemicals, glycyrrhetinic acid (GA) (Figure 2 B) ( (Noble et al., 2019)). In the case of GA, the reduction of relative calcium influx can be as large as 80%. It is important however to note that voltage (patch) clamp experiments do not show this action of SKT (Figure 2 C,D). We therefore hypothesise that the action by SKT on ICaL may be indirect, such as Calcium-induced inactivation of ICaL.
All simulations were run using COR or OpenCOR (version 0.6) (http://opencor.ws/) ( (Noble et al., 2019)). The original CellML file along with all the codes can be found in the following link in the PMR: https://models.physiomeproject.org/workspace/5f4 It is worth to mention that in order to be able to run sedml file the corresponding cellml file needs to be downloaded in the same folder. From (Noble et al., 2016).

Validation of the extended model, in which sarcolemmal calcium transporters are incorporated
For validation of the extended model, we used the data on the actions of SKT shown in Figures 1 and 2. We tested whether the extended model can reproduce the dynamics of actions of SKT when the L-type calcium channel is blocked by 80%, to correspond to the maximum block of calcium influx observed in experiments with the pure chemical GA (Figure 2, right).

3/16
The result shows a simple exponential decline in SR calcium, free intracellular calcium and the contraction ( Figure 3). The time constant of decline of the contraction was 22.2 minutes. This is very similar to the experimentally observed time constant of decline in contraction by 10 mg/ml SKT: around 21 minutes in Figure 1. The correspondence between the experimental and model results is very close. Therefore, the extended model succeeds in reproducing the dynamics of SKT when calcium influx is blocked by 80%. The simulation file Fig03.sedml contains the computational setting for running the model.

Sensitivity Analysis
It is important to determine how robust the extended model is. We carried out sensitivity analysis by increasing all the calcium regulatory network parameters by 50%, to test the sensitivity of the extended model to likely experimental variations in parameters. The results show that the extended model is robust. No changes in dynamics were greater than the change in parameter ( Figure 4). The correspondence between the experimental and the extended model results for the dynamics of contraction change following block of ICaL is therefore reliable. It is important to note that, apart from adding surface membrane transporters for calcium, we made no other changes to the original model.

Experiments in silico
Comparison of the extended and original models and reconstruction of muscle cramp.
The extended model can therefore be used to perform further experiments, to explore the possible reproduction of muscle cramp. First, we ran the extended model for 100 s at a stimulation frequency of 100 Hz, which produces a tetanic contraction that rapidly declines (fatigue), as illustrated by Figure 5. The simulation file Fig05.sedml contains the computational setting for running the model.  (Noble et al., 2019). Each calcium regulatory parameter was increased by 50%. Changes in the time constant were always less than 50%.
Figure 5. The extended model was run for 100 seconds at a stimulation frequency of 100 Hz to produce tetanic contractions leading rapidly to fatigue. The traces are: A -membrane potential, B -interstitial potassium, C -SR calcium, D -free intracellular calcium, E -L-type calcium current, and F -contraction (cross-bridge reaction). Note that the slow rise in contraction (actually cross-bridge reaction rate) is so large that the initial response appears very small. The rapid fatigue process is more clearly visible in the contraction panel of Figure 6. The results presented in this figure can be reproduced with the Fig05.sedml Prediction by modelling experiments: (a) accumulation of interstitial potassium, rising to nearly 10 mM. This magnitude of increase in muscle interstitial potassium has been recorded experimentally in human subjects during strong contractions ( (Green et al., 2000)) (b) gradual rise in SR calcium and calcium release (c) gradual rise in free intracellular calcium.

5/16
(e) large and gradual increase in contraction (rate of cross-bridge reaction).
Thus, the increase in contraction follows the increase in free intracellular calcium in time course, as shown in figure (c) and (e). Two questions then arise. First, how does this result differ from that of the original unextended model? To answer this question, we repeated the experiment using the same stimulation protocol as in Figure 5, but with ICaL completely blocked, corresponding to the fact that the original (Shorten et al., 2007) model has no ICaL. The result is shown in Figure  6. It is clear that the slow rise in free intracellular calcium and in contraction no longer occurs. The simulation file Fig06.sedml contains the computational setting for running the model. The simulation conditions are exactly same as the one in figure 5. We have also calculated the effect of reducing all surface calcium transporters to zero so that only the transporters in the original Shorten model are present. This is important since other calcium transporters, in addition to ICaL, were incorporated into the extended model. This is as expected since ICaL is responsible for nearly all calcium influx. The result is shown in Figure 7, which is identical to figure 6 As in Figure 6 there is no development of cramp. The simulation file Fig07.sedml contains the computational setting for running the model. The second question is whether the large slow rise in contraction could represent the cramp process? To answer that question, we performed the same experiment as in Figure 5, but switched the stimulus off after 100 seconds. To satisfy the definition of cramp, the contraction should not be entirely reversible when the stimuli cease. The result is shown in Figure 8. At 100 seconds, the magnitude of the contraction in the extended model rapidly drops, but the highly contracted state is not immediately reversible. The rapid drop is attributable to re-uptake by the SR (third trace). The contraction that remains, however, is still very large. It is around 0.5 units, which is 5 times larger than the peak tetanic contraction (0.1 units) before fatigue in Figure 6. This remaining strong contraction then slowly declines. Thus, the result in Figure   Lower frequencies of stimulation produced small degrees of calcium accumulation. By the definition of cramp we are using, frequencies of 25 Hz and lower did not produce the cramped state. We have therefore used 100 Hz as standard, which readily produced cramp and permits computational experiments on possible mechanisms of reversal. Figure 9 shows the result of a 10-minute run showing the free intracellular calcium and contraction on higher scales. The free intracellular calcium remains high throughout the run. The contraction remains above peak tetanic contraction for about 2 minutes, and remains high for about 8 minutes. By the usual criterion defining the condition of cramp, the results in Figures 8 and 9 show that the extended model automatically produces a model of cramp. The important discovery is that calcium influx is essential for producing cramp, as shown in Figures 6, 8 and 9. The simulation file Fig09.sedml contains the computational setting for running the model.

Reversal of the cramp condition
The comparison between Figures 5 and 6 shows that reducing calcium influx would be expected to delay or prevent cramp. But, how could cramp be reversed, once the cramped state has developed? To answer this question we performed the computational experiments shown in   Figures 8 and 10 shows that block of ICaL does not reverse the cramped state once it has developed.Because calcium influx is essential in the model for creating the cramp condition, it is reasonable to think that calcium efflux must be a necessary component of its reversal. We therefore performed computational experiments to explore the possible mechanisms that could produce sufficient calcium efflux for reversing the cramp state. We found that the extended model successfully reproduces the process of reversing cramp, provided that the accumulation of potassium is rapidly reversed. As we showed in the other paper, ( (Noble et al., 2019)), the extended model includes this process. Figure 10. This computational experiment uses the same protocol as Figure 7 except that the L-type calcium channel was blocked after stimuli were stopped (i.e. after 100 seconds). The traces are: A -membrane potential, B -interstitial potassium, C -SR calcium, D -free intracellular calcium, E -calcium current (ICaL), F -and contraction (cross-bridge reaction).
The results presented in this figure can be reproduced with the Fig10.sedml The result is shown in Figure 11. The protocol was the same as in Figures 8 and 10 up to the point where simulation was ceased (i.e. after 100 seconds). At that point: (a) the bulk extracellular potassium concentration was reduced to an assumed normal plasma level (4 mM).
(b) the time constants of exchange with bulk potassium and between tubular and interstitial potassium were reduced to 20% of their standard values to mimic rapid wash-out via increased blood flow.
The result in Figure 11 shows a rapid wash-out of the accumulated potassium, which is almost complete within 20 seconds. The cramp state is completely reversed within around 10 seconds. Wash-out of potassium is sufficient to reverse cramp condition. The simulation file Fig11.sedml contains the computational setting for running the model. Figure 12 shows superimposition of the contraction traces from the results illustrated in Figures  10 and 11. Without potassium wash-out, the cramp state persists for at least 2 minutes. With potassium wash-out, the cramp state disappears within 10 seconds. In order to investigate the interactions between the processes involved, we performed computational experiments on an expanded time scale, as shown in Figures 13 and 14. Figure 13 shows a more complete set of traces obtained on a larger time scale by calculating the changes occurring during the first 20 seconds of the simulation in Figure 10 following cessation of stimulation. The expanded time scale in Figure 13 enables us to investigate the interactions between the processes involved. Calcium efflux following potassium wash-out occurs through activation of sodium-calcium exchange.  The simulation file Fig13.sedml contains the computational setting for running the model. After 10 seconds after the simulation is stopped and values are changed, simulation results were cleared and the ending point was changed to 20000 milliseconds and all the stimulation was removed, finally simulation was run again for 20 seconds.

Superimposition of the cramped and relieved states
The results in Figure 14 shows that the fall of free intracellular calcium (fourth trace) is largely produced by activation of sodium-calcium exchange (bottom trace). The top two traces in Figure  14 show that the change in membrane potential in a hyperpolarizing direction closely follows the change in interstitial potassium, as expected since the resting potential is largely determined by the Nernst equation for potassium ions. Free intracellular calcium (fourth trace) initially falls very rapidly. The simulation file Fig14.sedml contains the computational setting for running the model. The simulation conditions are exactly same as the one in figure 13. Figure 13. This shows changes in A -membrane potential, B -interstitial potassium, C -SR calcium, D -free intracellular calcium, E -L-type calcium current , F -contraction and G -the sodium-calcium exchange current during the first 20 seconds after stimulation ceases in the computation illustrated in Figure 11. The results presented in this figure can be reproduced with the Fig13.sedml The changes are shown on a larger time scale in Figure 14, where we have also included the changes in the calcium pump current, which plays a relatively minor role in the model. The middle trace (contraction) in Figure 14 shows that there are two phases of fall in contraction in the model. The first, rapid, phase lasts about 200 milliseconds and is largely driven by the fall in free intracellular calcium produced by the high activity of the sodium-calcium exchange. This phase corresponds to the time taken to remove the calcium overload. Free intracellular calcium falls towards the usual maximum level during contraction, i.e. around 1 µM. There is then a slower phase of recovery driven by both the sodium-calcium exchange and the calcium pump.
The computations shown in Figures 5-13 were all performed with the time constant (TauK3) for exchange between interstitial [K] e and bulk [K] b set to a standard value, which was arbitrarily set to the value of the time constant for exchange between tubular and interstitial potassium in the original (Shorten et al., 2007) model, i.e.TauK2. Given the greater distances that must be involved in the exchange with [K] b , the time constant for exchange with bulk potassium must Figure 14. This shows changes in A -SR calcium, B -free intracellular calcium, C -Contraction, D -the sodium-calcium exchange current, E -the calcium pump current during the first 20 seconds after stimulation ceases in the computation illustrated in Figure. The results presented in this figure can be reproduced with the Fig14.sedml actually be larger than TauK2. Moreover, in a cramped muscle, circulation will be reduced. We therefore performed calculations with TauK3 set to 200 seconds to represent such a situation ( Figure 15). There are several important outcomes of this experiment: (a) a shorter period of stimulation (60 seconds) was needed to model the cramp process.
(b) accumulation of interstitial potassium was greater and so was the magnitude of the cramped state as a result.
(c) the contraction, after stimulation stopped, decayed more slowly and lasted much longer.
(d) when potassium wash-out was induced by reducing TauK3 to 4 seconds, the cramp state was relieved very rapidly.
It is noteworthy, however, that the increase in duration of the maintained contraction after stimulation ceases (about 3 fold) is not as long as the increase in time constant (about 10 fold). Another process must therefore be contributing to reducing interstitial potassium. The simulation files Fig15(01).sedml and Fig15(02).sedmlcontain the computational setting for running the model.
As noted in the other paper ( (Noble et al., 2019)), the initial steady state contraction in the extended model is very low compared to the original model. In that paper, we outlined how that deficiency might be rectified in future re-parameterisation of the extended model. In particular, we showed that increasing the calcium uptake rate, nuSR, into the SR can readily increase steady state contraction.  (02) We have therefore also run computations with the nuSR parameter increased to 15 and 35. Both levels of nuSR produce much larger steady state contractions, as shown in Figure 7 of the other paper ( (Noble et al., 2019)), but the development of the maintained (cramp) condition still occurs, as shown in Figure 16. The level of maintained contraction is much larger than in Figure 15, 1.25 compared to 0.35, but decays faster. The faster decay is attributable to a larger activation of the sodium-calcium exchange current by larger intracellular contraction. The simulation file Fig16.sedml contains the computational setting for running the model. The value of TauK3 was set to 200 seconds. The interplays between the various parameter adjustments that can restore larger steady state contractions in the model are complex and would require further investigation of the parameter spaces involved in future work.

Significance of sarcolemmal calcium fluxes in skeletal muscle
Sarcolemmal calcium influx is not necessary for individual twitch responses in skeletal muscle since the trigger for electromechanical coupling is entirely electrical. This property is what enabled Shorten et alShorten et al. (2007) to reproduce contraction and 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 tens of seconds the sarcolemmal calcium fluxes would be expected to have a significant effect on intracellular calcium concentrations and therefore on electromechanical coupling. We show in this paper that these small fluxes do gradually alter SR calcium and intracellular free calcium. The slow time course of those changes closely corresponds to the dynamics of change in contraction observed experimentally when a medication that reduces calcium influx is used. This close correspondence with experimental data partially validates the model.

Significance of calcium efflux, conditions for cramp, its prevention, and its reversal
One of the aims of our extension of the model is to reproduce cramp. We define cramp to be a maintained contraction that continues for a substantial period of time without stimulation and is therefore initially not reversible by nervous control.
The model satisfies this condition. Even though small on a beat-to-beat time scale, the increased calcium influx at high frequencies is predicted to be sufficient to produce calcium overload to generate a cramped state. The extended model is thus the first successful mathematical model of skeletal muscle cramp. It is useful as a model for the study of cramp induced by vigorous exercise. However, it remains to be seen whether rest cramps might also be reproduced by the model. A possible clue from our work is the facilitation and maintenance of cramp by increased interstitial potassium. If rest cramps are produced through poor circulation, leading to accumulation of interstitial potassium, then they might have a common mechanism. That question remains to be answered.
We show that the model also makes an important prediction on how cramp may be relieved rapidly. The prediction is that wash-out of interstitial potassium is necessary to relieve the cramp process. Effectively, high interstitial potassium maintains the cramped state. Relief depends on increased calcium efflux via the sodium-calcium exchange led by the wash-out of interstitial potassium. The model therefore indicates that a mechanism that relieves cramp must improve the circulation through vasodilation, so achieving wash-out of interstitial potassium. All these predictions are experimentally testable.
The result shown in Figure 14 implies that another process must also be contributing to re-

14/16
moving interstitial potassium. In the extended model, the only process that can do this is the sodium-potassium (Na + ,K + -ATPase) pump, which generates potassium influx in exchange for sodium efflux. This interpretation is strongly supported by Nielsen et al. ((Nielsen et al., 2004)) who showed the effect of high-intensity training on the kinetics of interstitial potassium changes in human subjects. Training increases the activity of the Na + ,K + -ATPase pumps and reduces the interstitial potassium accumulation during exercise.
Of course, it will require further investigations of the extended model to clarify the conditions under which this process develops and how it recovers. These model investigations would best be done in conjunction with human experiments on cramp following vigorous exercise, and to compare the model predictions with experimental recordings of biochemical parameters during vigorous exercise.

Further discussion from an evolutionary perspective
Our results also have implications for the evolution of skeletal muscle in having become so different from smooth and cardiac muscle in the role of calcium influx. Very low levels of calcium influx in skeletal muscle could be seen to be adaptive since high levels of calcium influx would generate the cramp mechanism too rapidly during maintained tetanic contractions. It is significant that we needed to use high-frequency stimulation to generate the cramp state.
Source code: a standard procedural language is used (e.g. MATLAB, Python, C) 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