Exhaustion of Skeletal Muscle Fibers Within Seconds: Incorporating Phosphate Kinetics Into a Hill-Type Model

Initiated by neural impulses and subsequent calcium release, skeletal muscle fibers contract (actively generate force) as a result of repetitive power strokes of acto-myosin cross-bridges. The energy required for performing these cross-bridge cycles is provided by the hydrolysis of adenosine triphosphate (ATP). The reaction products, adenosine diphosphate (ADP) and inorganic phosphate (P_i), are then used—among other reactants, such as creatine phosphate—to refuel the ATP energy storage. However, similar to yeasts that perish at the hands of their own waste, the hydrolysis reaction products diminish the chemical potential of ATP and thus inhibit the muscle's force generation as their concentration rises. We suggest to use the term “exhaustion” for force reduction (fatigue) that is caused by combined P_i and ADP accumulation along with a possible reduction in ATP concentration. On the basis of bio-chemical kinetics, we present a model of muscle fiber exhaustion based on hydrolytic ATP-ADP-P_i dynamics, which are assumed to be length- and calcium activity-dependent. Written in terms of differential-algebraic equations, the new sub-model allows to enhance existing Hill-type excitation-contraction models in a straightforward way. Measured time courses of force decay during isometric contractions of rabbit M. gastrocnemius and M. plantaris were employed for model verification, with the finding that our suggested model enhancement proved eminently promising. We discuss implications of our model approach for enhancing muscle models in general, as well as a few aspects regarding the significance of phosphate kinetics as one contributor to muscle fatigue.

1. Introduction

Muscles performing mechanical work become exhausted, that is, they fail to maintain high force levels for a longer time period. What seems like the most obvious statement for everyone conducting physical exercises has not found its way into large parts of biomechanical modeling. Several state-of-the-art muscle models do not comprise the physiologically well-observed force decay over time, especially in fast-twitch fibers under high neural stimulation (Brown and Loeb, 2000; Rode et al., 2009; Blümel et al., 2012; Millard et al., 2013; Haeufle et al., 2014b; Schappacher-Tilp et al., 2015; Mörl et al., 2016). Moreover, literature terminology is vaguely summarizing any possible force-decay mechanism under the umbrella of fatigue, which is discussed throughout disciplines, such as biology (Enoka and Stuart, 1992; Fitts, 1994), biomechanics (Bergström and Hultman, 1985; MacIntosh et al., 2012), engineering (Böl, 2009), medicine (Chen et al., 1999; Cardozo et al., 2011), physiology (Lindström et al., 1977; Allen and Westerblad, 2001; Westerblad and Allen, 2002; Allen et al., 2008), and sports science (Komi, 2000; Brown et al., 2017). One of the first books addressing (muscle as well as whole body) fatigue was published at the beginning of the 20th century (Mosso, 1904) and a multitude of research has followed since, see Gandevia (2001) for a thorough review.

Commonly, muscle fatigue, i.e., the decline of the generable force level over time, is differentiated between central fatigue, i.e., the inability of the neural network to provide sufficient stimulation, and peripheral fatigue, i.e., the inability of the muscle cells to provide energy through metabolic activities (cf. Bigland and Lippold, 1954; Vøllestad, 1997). The main part (ca. 80%) of force decay is hereby associated with the latter, metabolic component (Kent-Braun, 1999). Due to the broad scientific interest, many metabolic factors of force decay have been identified, for example adenosine triphosphate (ATP) consumption (Sieck et al., 2013), calcium precipitation (Baylor and Hollingworth, 1998; Allen et al., 2008), glycolysis (Abbiss and Laursen, 2005), intracellular acidosis (Gandevia et al., 1995, Ch. 3), lactic acid (Tesch et al., 1978; Westerblad et al., 2002), pH value (Westerblad and Allen, 2002; Cooke, 2007), potassium loss or accumulation (Bickham, 2003; Allen et al., 2008) (Gandevia et al., 1995, Ch. 5), and reactive oxygen species (Debold et al., 2004; Cooke, 2007; Allen et al., 2008; Glass, 2017). A summarizing flowchart can be found in Abbiss and Laursen (2005, Figure 9). Although the role of some factors is controversial, e.g., the role of lactic acid (Tesch et al., 1978; Westerblad et al., 2002) or protons in general (Debold et al., 2016), it is agreed upon that the accumulation of inorganic phosphates (P_i) accounts for the main factor in peripheral fatigue (Hibberd et al., 1985; Nosek et al., 1987; Kowaltowski et al., 1996; Allen and Westerblad, 2001; Westerblad et al., 2002; Debold et al., 2004). Despite vast physiological findings, biomechanical models, including force decay on the time scale of seconds, focus on rather phenomenological descriptions of typical force decay patterns (Liu et al., 2002; El ahrache et al., 2006; Xia and Frey-Law, 2008; Ma et al., 2009, 2012; Frey-Law et al., 2012; Rashedi and Nussbaum, 2015b), or try to cover the entirety of physiological processes (Shorten et al., 2007) which is computationally expensive.

Here, we develop a mathematically straightforward, yet physiology-based model that is able to explain the majority of early force decay while being computationally manageable. We introduce a model for hydrolytic ATP-ADP-P_i dynamics (or short phosphate dynamics) (Allen and Orchard, 1987) and the relative chemical potential ${\tilde{\mu }}_{\text{ATP}}$ of ATP as the corresponding state variable that is diminished by both ATP depletion as well as P_i accumulation. As validation, isometric contraction data at different muscle lengths from rabbit M. gastrocnemius and M. plantaris were used to reproduce typical force-time patterns. Our aim is to provide a basic enhancement of the well-known realm of Hill-type muscle models by describing bio-chemical kinetics that can be altered or extended, depending on the experimental resolution at hand.

2. Muscle Data

Experiments on female New Zealand white rabbits (Oryctolagus cuniculus, n = 2) with an average weight of 3.06 ± 0.01 kg (mean±SD) were approved according to section 8 of the German animal protection law (Tierschutzgesetz, BGBl. I 1972, 1277). The animals were anesthetized with Bupivacain (Jenapharm, 1 ml, 0.5%, epidural) after short-term sedation with sodium pentobarbital (Nembutal, 80 mg/kg body weight) and all efforts were made to minimize suffering. Anesthesia, preparation as well as experimental setup have been previously described (Böl et al., 2013, 2015; Siebert et al., 2015). In brief, the muscle to be examined was freed from its surrounding tissues and the rabbit was fixed with bilateral bone pins to a stereotaxic frame. The distal tendon of the muscle was attached horizontally to a muscle lever system (Aurora Scientific 310B-LR). Because the muscle performance depends on temperature, the animal was heated during the complete experimental procedure using a heating pad (Harvard Apparatus, 39.0 ± 0.4°C, mean±SD) and the surface of the isolated muscle was frequently sprinkled with heated (39°C) physiological saline solution.

For this study, we particularly investigated the performance of M. gastrocnemius (GAS) and M. plantaris (PLA). At the beginning of all measurements, the corresponding initial muscle lengths were determined in situ with a micrometer at an ankle and knee joint angle of 90°. Then a series of isometric experiments was performed at different muscle lengths (GAS: 114–120 and 128–132 mm; PLA: 114–124 mm; length increments of 2 mm) comprising the force-length relation from the beginning of its ascending limb (about 10% maximum isometric force, F_max) to its plateau region (100% F_max). To avoid muscle damage, the muscles were lengthened until passive forces reached about 10 and 20% of the maximum force for PLA and GAS, respectively. Thus, for GAS, there are two additional isometric measurements at the beginning of the descending limb of the force length relation (130 and 132 mm muscle length). Muscles were stimulated (Aurora Scientific 701C) with 100μs square wave pulses at 130 Hz (supramaximal tetanic muscle stimulation) via the tibial nerve using a bipolar gold electrode for 700 ms. One additional measurement for the GAS at 126 mm was conducted with 1,400 ms of stimulation.

3. Incorporating Bio-chemical Kinetics Into a Macroscopic Muscle Model

In this section, we work out an enhancement of an existing Hill-type muscle model (see Appendix A), which describes the dynamics of the activation and contraction of muscle fiber material. The model is enhanced by an additional process (phosphate dynamics), which is formulated in terms of a differential-algebraic equation.

Beforehand, in literature, we identified five distinguishable approaches of modeling short-term force decay patterns in skeletal muscle as a consequence of ongoing stimulation. Arranged with respect to increasing physiological verisimilitude these are

(1) Functional effects (multi-parametric, control theoretical): introducing a muscle compartment model, whose states can be controlled to match desired data trajectories (Hawkins and Hull, 1993; Freund and Takala, 2001; Liu et al., 2002; Xia and Frey-Law, 2008; Böl, 2009; Frey-Law et al., 2012; Sih et al., 2012). For computational simplification, Wiener-Hammerstein models (Wiener, 1942; Narendra and Gallman, 1966; Wills et al., 2013) are also used (Cai et al., 2009).

(2) Isolated functional effects (low-parametric): search for a functional dependency of force decay on calcium ions (Dorgan and O'Malley, 1998), calcium-troponin complexes (Ding et al., 2000; Marion et al., 2010), pH value (Giat et al., 1993), or glycolytic flux (Callahan et al., 2016) partly involving purely descriptive parameters bearing no relation to physiology.

(3) Phenomenological details: introducing a theoretical construct, named for example fatigue index (Ma et al., 2009, 2011), fitness level (Tang et al., 2005), fatigue rate (Ma et al., 2012), co-contraction factor (Seth et al., 2016) or similar (Deeb et al., 1992; James and Green, 2012), which describes the empirically observed load-endurance relation per (linear) ODE. Reviews of applying this method can be found in (El ahrache et al., 2006; Ma et al., 2009; Rashedi and Nussbaum, 2015a,b).

(4) Phenomenological chains: mathematically describing macroscopic processes based on measured aerobic and anaerobic power output (Eriksson et al., 2016).

(5) Resolved physiological chains: mechanistically formulating interactions of bio-chemical processes with predictive power; potentially resulting in dozens of (partial) differential equations and more than 100 parameters (Shorten et al., 2007).

We aim at predicting the force decay of a biomechanical muscle computer model as response to a given muscle stimulus (forward dynamic simulation). Opposing, in approach (1) the force level is considered given and the stimulus is to be estimated in order to match the observed force decay (inverse dynamic simulation). Thus, here approach (1) can be ruled out. As well can approach (4), because it does not incorporate inner-muscular properties. An advantage of approaches (2) and (3) is the computational and epistemological simplicity, whereas a disadvantage is the missing physiological interpretability of the model parameters. In contrast, the physiological backbone of approach (5) is beyond dispute, however, so is its computational complexity. Our aim is to combine the strengths of approaches (2), (3), and (5), simultaneously neglecting their weaknesses, by developing a computationally cheap model based on physiological knowledge (cf. also section 5.1). As mentioned before as well as in Shorten et al. (2007), the accumulation of P_i as a consequence of ATP turnover accounts for the main reason for short-term force decay. Therefore, a single, linear ODE is added to the existing model (see Equation 3), describing the change in [ATP] and [P_i] within the sarcoplasmatic reticulum and its interplay with cross-bridge mechanics, based on known reaction kinetics, for a detailed description see Appendix B.

Following the iconic paper of Lymn and Taylor (1971), ATP binds to the attached myosin head to allow its release from the actin filament. The energy release in the hydrolysis is then used to re-configure the myosin head for re-attachment in order to undergo another power stroke. Although many physiological details remain uncertain, the core element of the cross-bridge cycle can be represented by a simple pseudo-first order equilibrium

$\begin{array}{ll}\text{ATP}\underset{k_{\text{con}}:\text{ ATP Synthase}-{\text{H}}_2\text{O}}{\overset{k_{\text{hyd}}:\text{ ATPase}+{\text{H}}_2\text{O}}{\rightleftharpoons }}\text{ADP}+{\text{P}}_i,& (1)\end{array}$

namely ATP hydrolysis/condensation. Herein, k_con and k_hyd denote the rates of condensation (catalyzed by mitochondrial ATP synthase) and hydrolysis (catalyzed by ATPase within the myosin head), respectively, with the corresponding equilibrium constant

$\begin{array}{ll}{K}_{\text{ATP}}=\frac{k_{\text{con}}}{k_{\text{hyd}}}=\frac{[\text{ATP}]·c_0}{[\text{ADP}]·[{\text{P}}_i]},& (2)\end{array}$

where c₀ = 1M represents the standard concentration serving as a normalization factor (Allen and Orchard, 1987; Hancock et al., 2005; Cooke, 2007). Following Equation (1), the time evolution of [ATP] or [P_i] can be expressed in terms of the first order ODE

$\begin{array}{ll}\frac{\text{d}}{\text{d}t}[{\text{P}}_i]=k_{\text{hyd}}·[\text{ATP}]·c_0-k_{\text{con}}·[\text{ADP}]·[{\text{P}}_i].& (3)\end{array}$

Note that although this dynamic is motivated by physiological considerations, it constitutes for a vastly oversimplified description of the bio-chemical processes underlying the full cross-bridge cycle. However, Equation (3) can be enhanced if explicit measurements (e.g., concentrations or time constants) of the involved reactants were available, see section 5.2 and in particular Figure 5.

For a start, in our model, K_ATP and in particular k_hyd, are assumed to be dependent on the current amount of calcium-bound troponin C terminals $\tilde{q}$ as well as relative fiber length ${\tilde{\ell }}_{\text{CE}}$, cf. Appendix A. The condensation rate k_con is assumed to be a constant. The dependence of k_hyd on $\tilde{q}$ is obvious, because ATP hydrolysis within the myosin head does not occur if there are no available active sites on actin. It has been shown that, within resting muscle fibers, there is ATP consumption, but far less than in activated ones (Hilber et al., 2001, Figure 2A). This resting utilization corresponds to our model parameter ${\tilde{q}}_{min}$ representing minimum troponin-activity. Hence, the rate constant k_hyd is assumed to increase along with $\tilde{q}$, i.e., to shift the reaction Equation (1) in favor of the hydrolysis products. Experiments further show that both ATP consumption rate (Aljure and Borrero, 1968; Doud and Walsh, 1995; MacNaughton and MacIntosh, 2007; Fenwick et al., 2016) and P_i release (Bickham et al., 2011, Figure 3a) are inversely proportional to fiber length (see also section 5.5). Summarizing, the following ansatz is obtained:

$\begin{array}{ll}{K}_{\text{ATP}}=K_{\text{ATP}}(\tilde{q},{\tilde{\ell }}_{\text{CE}})=\frac{k_{\text{con}}}{k_{\text{hyd}}(\tilde{q},{\tilde{\ell }}_{\text{CE}})}=\frac{k_{\text{con}}}{{\widehat{k}}_{\text{hyd}}·\tilde{q}/{\tilde{\ell }}_{\text{CE}}},& (4)\end{array}$

with the parameter ${\widehat{k}}_{\text{hyd}}$ representing the hydrolysis rate at full activity and optimal fiber length. Note that the alterations of the rate constants with respect to external conditions, such as fiber type, pH value, or temperature, are not (yet) incorporated (cf. section 5.3). It is further assumed that in a sufficiently activated fiber there is enough calcium to simultaneously enable the ATPase within the myosin head (De La Cruz and Ostap, 2009; MacIntosh et al., 2012; Glass, 2017). Note that magnesium, which serves a co-factor, is not yet incorporated within our model, and assumed to be sufficiently available.

Ultimately, two mechanisms of phosphate dynamics are incorporated

(a) ATP is to some extent refueled by phosphate storages, such as Creatinephosphate and ADP, cf. Equations (22)–(24) as well as Carlson and Wilkie (1974), Hilber et al. (2001), and

(b) a rising level of P_i (and simultaneously of ADP) deteriorates the relative chemical potential of ATP, denoted ${\tilde{\mu }}_{\text{ATP}}$, which is defined as

$\begin{array}{ll}{\tilde{\mu }}_{\text{ATP}}:=\frac{-\Delta G_{\text{ATP}}^{\circ }+R·T·ln\left(\frac{[\text{ATP}]·c_0}{[\text{ADP}]·[{\text{P}}_{\text{i}}]}\right)}{{\mu }_{\text{ATP,max}}},& (5)\end{array}$

where R denotes the universal gas constant, T the temperature, and μ_{ATP, max} ≈ 60 kJ/mol the maximum chemical potential of ATP in a resting fiber (Barclay, 2015) (see also Equations 37 and 38).

Incorporating the new state ${\tilde{\mu }}_{\text{ATP}}$ in our model, the term muscle activity (ã, cf. Appendix A) from now on refers to the relative amount of force-producing cross-bridges formed at the available active sites, i.e.,

$\begin{array}{ll}\tilde{a}:=\tilde{q}·{\tilde{\mu }}_{\text{ATP}}.& (6)\end{array}$

Compare Appendix A.0.3 and Equation 21 for a quick assessment of the changes that arise in the underlying model from Günther et al. (2007).

Certainly, our introduced model enhancement is far from capturing every physiological cause of force decay in skeletal muscle. Yet, given that the full mechanisms of the remaining factors are at best partially known (MacIntosh and Rassier, 2002; Allen et al., 2008), it may arguably serve as a physiological and partly mechanistic basis for further enhancements. As the results in section 4 show, the addition of a single ODE containing only one new parameter ${\widehat{k}}_{\text{hyd}}$ is already sufficient to model the short-term force decay in isometric contraction experiments. Table 1 summarizes the parameters necessary to describe the full phosphate kinetics as explicitly explained in Appendix B.

TABLE 1

Table 1. Overview of the nine new model parameters, necessary to describe phosphate dynamics (cf. Appendix B).

4. Results

In total, n_GAS = 7 and n_PLA = 6 isometric force-time datasets from Siebert et al. (2015) for GAS and PLA, respectively, were compared with our model's output when applying the same stimulation protocol (control function) as within the experiments. Following Figure 1, each experiment had a runtime of 1.2 s, with a sampling rate of 1 kHz. The control function is zero, i.e., no electrical stimulation, for the first 0.1 s and is set to one, i.e., full stimulation, for further 0.7 s. Within the last 0.4 s, the muscle again experiences zero stimulation and finds its pre-stimulation equilibrium. It had been priorly shown that such isometric experiments are generally utilizable to estimate dynamic parameters (Rockenfeller and Günther, 2016), particularly within the force rise and fall phases at the beginning and the end of the stimulation. Valid estimations for static parameters, such as slack lengths and maximum force, are possible at the equilibrium states for zero and full stimulation, respectively.

FIGURE 1

Figure 1. Force-time courses of GAS (left column) and PLA (right column) at various muscle lengths (see section 2). Models were fit to data, where different lengths are indicated by different colors as labeled in the top row. Parameter values are summarized in Table 2. (A,B) Classical model (GAS⁻, PLA⁻) without phosphate dynamics. (C,D) New model (GAS⁺, PLA⁺) including phosphate dynamics. (E,F) New model with switched off ATP hydrolysis (GAS^{+, 0}, PLA^{+, 0}, i.e., ${\widehat{k}}_{\text{hyd}}=0$).

In order to assess the effect of the new state introduced in Equations (3) and (37), as well as to increase validity, a comparative parameter estimation was conducted; once with the basic model as presented in Appendix A and once with the enhanced model. The formalization of the underlying optimization procedure can be found in Appendix C. Let GAS⁻ and PLA⁻ denote the models without phosphate dynamics, and let GAS⁺ and PLA⁺ indicate the inclusion. Further, we show the effect of switching off phosphate dynamics (${\widehat{k}}_{\text{hyd}}=0$) in the inclusive model for which the nomenclature GAS^{+, 0} and PLA^{+, 0} was used. Table 2 lists the obtained parameter values as well as the pre-set boundaries (see next section 4.1). Figure 1 shows the associated force-time model outputs. These force-time results, obtained by least-square fits, are compared with respect to residues for both L₁ (absolute) and L₂ (squared) distances for the purpose of showing consistency throughout the model variants.

TABLE 2

Table 2. Bounds and optimized parameter values resulting from the estimation process, described in Appendix C.

To further support the validity of the exhaustion model, we compared the model prediction with one additionally available experimental data set, namely an isometric contraction of GAS with longer stimulation time (1.4 s) at medium length (126 mm) (see Figure 2). In agreement with experiments, the model simulation yielded a decay in muscle force of 23% over the stimulation period. Thus, the model is able to predict isometric experiments with prolonged stimulation times. However, this statement only holds true for a single additional isometric contraction. Further long-stimulation experiments, also with PLA, should be conducted in order to strengthen this claim. As solely determined by isometric contractions, the herein obtained model parameters should in any case be cautiously considered when aiming at dynamic simulations with high shortening velocities.

FIGURE 2

Figure 2. Force-time course of an additional isometric experiment on GAS (dots) as well as the model output of GAS⁺ (gray line), utilizing parameter values from Table 2. Stimulation duration is 1.4 s, GAS length ℓ_GAS = 126mm. Inlay shows the concentration-time courses of the most relevant bio-chemical reactants (blue: ATP; red: P_i; green PCR) as well as the relative chemical potential of ATP (black: ${\tilde{\mu }}_{\text{ATP}}$). For more details see Figure A1 in Appendix B.

4.1. Parameter Estimation

4.1.1. Setting the Parameter Bounds

As presented in Appendix C, parameter estimation/optimization was conducted in a least-squares sense. It has been further possible to restrict the search area for the algorithm, i.e., pre-define a hyperrectangle in the parameter space by lower and upper bounds, see third and sixth column of Table 2. Most parameters settle well in between these pre-set bounds, which were chosen based on physiological plausible values and prior model knowledge (Mörl et al., 2012; Siebert et al., 2015; Rockenfeller and Günther, 2016). As far as possible, parameter bounds were given similar for GAS and PLA. Naturally, deviations in parameter values for these different muscles occurred in (slack) lengths as well as forces. Further, a challenge occurred due to the redundancy in serial and parallel elasticities. As is apparent from Equation (9), serial and parallel elastic element are assumed to operate in series. Further, MTU length is held constant (static) in isometric experiments. In consequence, the individual contributions of parallel and elastic elasticities of the muscle in dynamic situations can not be clearly resolved. Hence, restrictions in either parameter set might influence the convergence property of the other and the herein presented bounds have to be considered with caution.

4.1.2. Comparison Between Model⁻ and Model⁺

The most prominent effect of including phosphate dynamics on parameter values is an about 10% increase of the parameter value of maximum isometric force (F_max) as calculated by the optimization algorithm. The explanation of an increasing F_max is straight forward: The model⁺ has to compensate the force decay induced by phosphate accumulation. This effect becomes particularly visible in the terminally investigated model^{+, 0}, where ATP hydrolysis and thus phosphate accumulation is artificially set to zero (see Figures 1E,F and next paragraph). Accordingly, the maximum velocity decreases ~40% in the model⁺, mostly due to a decreasing B_{rel, 0}. Finally, serial damping is also less in the model⁺.

For the rest of the parameters, there is no systematic change detectable. Optimal fiber length increases from GAS⁻ to GAS⁺ but decreases from PLA⁻ to PLA⁺. The opposite holds true for the serial elastic element (SEE) parameters ΔU_{SEE, nll} and ΔU_{SEE, l}. Activation dynamics parameters remain almost constant.

4.1.3. Comparison Between Model⁺ and Model^{+, 0}

When switching off ATP hydrolysis in the model⁺, thus inhibiting phosphate accumulation, model force continues to rise well above measured force levels, particularly at longer MTU length (cf. Figures 1E,F). This observation is in accordance with experiments (Phillips et al., 1993), where a decrease of inorganic phosphate was associated with an increase in force production (Westerblad and Allen, 2002). Hence, the value of F_max in the model⁺ now contains additional information about a theoretically achievable maximum force, if inorganic phosphate was pumped out of the muscle cells instead of accumulated.

4.1.4. Comparison Between GAS and PLA Models

The optimized parameters are in good agreement with prior experiments (Siebert et al., 2015) (see also section 5.4). In only two (GAS), respectively one (PLA) experiment the CE reaches lengths above ℓ_{CE, opt} and forces close to F_max. Hence, the shape of the descending branch of the force-length relationship as well as the results for the non-linear to linear transition of the SEE have generally to be taken with caution (cf. Figure 4).

In contrast to GAS⁻, the PLA⁻ model shows a more or less pronounced force overshoot after being fully stimulated, see particularly Figure 1B at medium lengths between 0.1 and 0.2 s. While trying to compensate for the apparent force decay, the optimizer found an interesting, non-trivial parameter setup: serial damping was adjusted to become very strong (D_SDE = 10) and completely force independent (R_SDE = 1), cf. (Günther et al., 2007, Equation 23). In return, the curvature of the eccentric part of the force-velocity relation was substantially increased (S_e = 4.5) in order to compensate a slow force decay at the end of the stimulation.

Maximum shortening velocity (B_{rel, 0}/A_{rel, 0}) is twice as large in GAS⁻ compared to PLA⁻, i.e., 80–44 ℓ_{CE, opt}/s, and 50% larger for GAS⁺ compared to PLA⁺, i.e., 41–27 ℓ_{CE, opt}/s. Mainly this difference is explained by A_{rel, 0} being 3.5 times larger for GAS than for PLA (see Table 2).

The force-length relation of GAS shows a broad plateau-like region with steep ascending and descending limbs, whereas PLA shows a shallow ascending limb (see Figure 4). We already mentioned the lack of trustfulness of the descending limb's shape. However, the difference in the force-length relation might additionally be influenced by pennation of the fibers. As the fiber shortens, its pennation angle increases (Drazan et al., 2019). At higher forces and lower velocities, the force-length characteristic of the fiber is then transformed (geared) to an altered force-length characteristic of the whole muscle (Azizi et al., 2008).

Finally, PLA (${\widehat{k}}_{\text{hyd}}=2.27$) seems to metabolize ATP, and thus accumulate phosphate, more rapidly than GAS (${\widehat{k}}_{\text{hyd}}=1.45$), which is in good agreement with the fiber type composition of PLA (> 90% fast twitch fibers) and GAS (> 75% fast twitch fibers) (Wang and Kernell, 2001; Siebert et al., 2015).

4.1.5. Figures and Residues

Figure 1 shows the force-time measurements compared to the model evaluations for GAS (left column) and PLA (right column). Within these columns, the top row shows the corresponding best fit model⁻, the middle row the newly proposed model⁺ and the bottom row the effect of an a-posteriori neglect of exhaustion, i.e., model^{+, 0}. The residues in the optimized least-square sense (L₂) were for comparability scaled to 1000 data points and account for 61.43, 56.17, and 166.79 N for GAS⁻, GAS⁺, and GAS^{+, 0}, as well as 32.19, 28.58, and 70.26 N for PLA⁻, PLA⁺, and PLA^{+, 0}, respectively. Hence, fits for the model⁺ resulted in around 10% less error than for model⁻. Fits of GAS and PLA were comparably good with respect to the difference in maximum force, that is F_max was estimated twice as large for GAS than for PLA, and so were the residues. In order to give an alternative error estimate, independent of the objective function, the absolute (L₁) deviation of model and measurement scaled to a single data point was calculated. These values are 1.57, 1.42, and 4.42 N for GAS⁻, GAS⁺, and GAS^{+, 0} as well as 0.84, 0.71, and 1.19 N for PLA⁻, PLA⁺, and PLA^{+, 0}, respectively, therefore revealing the exact same qualitative behavior.

Summarizing, models⁺ yielded a superior fit, which additionally captured the force decay characteristic predicted in Equation (8), particularly for GAS⁺, i.e., small decay rate at short CE lengths increasing with length but becoming smaller again at CE lengths above ℓ_{CE, opt}, see Figures 1C, 4. In Figure 1B the aforementioned effect of the strong damper in terms of an early force overshoot can be observed, especially for short lengths.

4.2. Parameter Sensitivity

Figure 3 shows the time courses of the sensitivity for all model parameters for two exemplary cases, one for a long GAS muscle and one for a short PLA. The (local) sensitivity values (see again Appendix C) indicate how small changes in the parameter value would affect the model output at the corresponding time instances. The absolute values were for comparability truncated at a value of one, which can be interpreted as for example a 10% change of the parameter value results in a 10% change in model output. Altogether, Figure 3 may help the reader to estimate the influence of parameters across time and magnitude as well as to acknowledge the setting of bounds in Table 2.

FIGURE 3

Figure 3. Example sensitivity-time matrices for each model parameter at a long length GAS (A, ℓ_GAS = 132mm, cf. dark blue curve in Figure 1A) and a short length PLA (B, ℓ_PLA = 116mm, cf. yellow curve in Figure 1B). Large absolute sensitivities indicate the local importance of parameter value accuracy. For further explanation see text. For scaling purpose, sensitivity values had been truncated at absolute values >1. For example for ℓ_{SEE, 0}, sensitivities go up to absolute values in the magnitude of 10⁴.

The model is most sensitive to reference and slack lengths of CE and SEE, respectively. On the one hand, they are properly determinable, but on the other hand they have to be known with high accuracy (cf. Table 2). For the longer muscle (Figure 3A), SEE, PEE and descending branch parameters gain importance, since the CE operates at lengths above ℓ_{CE, opt}. Likewise, switches in sign of the sensitivities can be observed for all SEE and most PEE parameters, indicating the different influences in passive and active muscle states. For the shorter muscle (Figure 3B), parameters describing the ascending branch are of importance, as would have been expected. The influence of eccentric parameters is visible for the regions of force decay, be it due to phosphate accumulation or end of stimulation. The influence of Hill parameters on the model output is clearly visible after start and end of stimulation. Parameters for activation dynamics show similar behavior across MTU lengths; ${\tilde{q}}_{\text{min}}$ co-determines the passive force, whereas the remaining parameters begin to influence the model output right after the start and until the end of the stimulation.

4.3. Empirical and Theoretical Force Decay

As a direct consequence of Equations (6) and (21), the MTU force at the isometric steady-state ($\frac{\text{d}}{\text{d}t}\tilde{q}=\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0$) can be written as

$\begin{array}{ll}{F}_{\text{MTU}}{\mid }_{\frac{\text{d}}{\text{d}t}\tilde{q}=0,\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0}=\tilde{q}·{\tilde{\mu }}_{\text{ATP}}·F_{max}·{\tilde{F}}_{\text{isom}}+F_{\text{PEE}},& (7)\end{array}$

which at full troponin-activity ($\tilde{q}=1$) results in a force decay rate of

$\begin{array}{ll}\frac{\text{d}}{\text{d}t}{F}_{\text{MTU}}{\mid }_{\tilde{q}=1,\frac{\text{d}}{\text{d}t}\tilde{q}=0,\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0}=\frac{\text{d}}{\text{d}t}{\tilde{\mu }}_{\text{ATP}}·F_{max}·{\tilde{F}}_{\text{isom}}~{\tilde{F}}_{\text{isom}}/{\ell }_{\text{CE}},& (8)\end{array}$

see also Equations (37)–(39) for the latter proportionality. Indeed, the connection $\frac{\text{d}}{\text{d}t}{F}_{\text{MTU}}~{\tilde{F}}_{\text{isom}}/{\ell }_{\text{CE}}$ has, to our knowledge, not been formulated yet in the literature, but can be found in experimental data from rabbit psoas muscle (Hilber et al., 2001, Figure 2B) as well as in electromyography studies on center frequencies (Doud and Walsh, 1995, Figure 4). The former source omitted a functional dependency, whereas the latter gave a decreasing linear fit, although the shape of the altered sarcomere force-length relation was prominently shaped as described. Figure 4 shows the force rates of our data, compared to the theoretically predicted relation.

FIGURE 4

Figure 4. (A) Experimental force rates (time derivatives) of GAS at different MTU lengths, normalized to maximum force. Color encoding is the same as in Figure 1. Data were smoothed by a moving average filter with 40 ms width for clarity of depiction, but not for calculations. Inlay shows the mean force rates (solid lines) in the 0.1 s interval between 0.58 and 0.68 s, where near steady-state conditions were assumed. (B) Experimental force decay rates vs. CE length at different MTU lengths (colored circles) as extracted from (A) in comparison to theoretical modified CE force-length relation (Equation 8, black line). An additional data point (diamond) was taken from another GAS experiment (see Figure 2). The best-fit parameter κ_GAS = 0.163 s⁻¹ was found by an optimization routine, all other parameters were taken from Table 2. The unmodified and scaled force-length relation [${\kappa }_{\text{GAS}}·{\tilde{F}}_{\text{isom}}({\ell }_{\text{CE}})$, gray line] is given for reference. In (C,D) the procedure is mirrored for PLA, where κ_PLA = 0.165 s⁻¹.

5. Discussion

5.1. Adding Submodels of Disregarded Physiological Processes

To predict the decay in muscle force during isometric contractions, we have added a (third) ODE (Equation 3), which describes the dynamics of [ATP], to an existing Hill-type muscle model which consisted of already two differential and several algebraic equations (see Appendix A). A. V. Hill's hyperbolic force-velocity relation (Hill, 1938), empirically found for muscle fibers, is in the core of Hill-type models. Therefore, such models are of empirical, macroscopic, and consequently reduced character. They often show deficiencies in their capabilities of reproducing the wealth of physiological and experimental conditions. For compensating one such deficiency, we have herein followed the methodological path of step-wise enhancing a Hill-type model. There are strong indications that the Hill relation, which has been inferred from the combination of mechanical and thermodynamic measurements, is founded in structural properties of muscle fibers (Günther and Schmitt, 2010; Rosenfeld, 2012; Seow, 2013; Rosenfeld and Günther, 2014; Günther et al., 2018), which also means that the Hill relation is not restricted to steady-state contractions (Piazzesi et al., 2002; Lemaire et al., 2016; Rockenfeller and Günther, 2016). Furthermore, because the Hill relation originates in both mechanics and thermodynamics, it may already well represent basic properties of active muscle tissue during dynamic interactions with other body tissues. We therefore conclude that the methodological path chosen has a sound basis.

In this study, we have now enhanced the activation dynamics part of our initial model. In its initial formulation the model considered activation dynamics introduced by Hatze (1977, 1978, 1981) in a recently simplified variant (Rockenfeller and Günther, 2018). Hatze had thoroughly described the effect of a neural impulse on calcium dynamics up to the binding to troponin and subsequent clearance of tropomyosin from the actin helix (cf. processes P1 to P5 in Rockenfeller and Günther, 2016, Appendix A). The symbol $\tilde{q}$ for the troponin activity thus describes the relative amount of available (cleared) active sites on actin at a given filament overlap. However, he assumed that each possible cross-bridge would instantaneously form and generate force in the wake of troponin activation, thereby ignoring any further mechanisms and processes delaying or interfering with a cross-bridge's force generation. For example, as already motivated and sketched in Figure 5 in more detail, ATP hydrolysis reaction kinetics and phosphate accumulation are well-known mechanisms to have a bearing on cross-bridge dynamics. Therefore, we have introduced the ATP's chemical potential ${\tilde{\mu }}_{\text{ATP}}$ which is a transform of the additional state variable [ATP] that represents any force-reducing effects on single cross-bridges. This new state can likewise be interpreted as an attachment-to-detachment ratio in the sense of Huxley (1974).

FIGURE 5

Figure 5. Schematic cross-bridge (Lymn-Taylor) cycle (1-2-3-4-5-6-7-1). New abbreviations are used as follows: Troponin C (TnC), actin (A), open and closed myosin switch-2-region (M^O and M^C), light chain domain (LCD), catalytic domain (CD). For their explanations, see text. Although most processes are easily reversible, for the sake of clarity only one direction is depicted.

It does not seem expedient to aim here for a higher (structural and parametric) resolution of the underlying processes within our model of activation dynamics, given and solely based on our experimental set-up (isometric whole muscle contractions). Yet, the chosen methodological path of step-wise model enhancement by adding structure-based equations is perfectly designed to do so, namely, to factor in, on the basis of experiments, processes like receptor-ligand binding processes and state transitions (Bagshaw and Trentham, 1973; Trentham et al., 1976; Stein et al., 1981; Grigorenko et al., 2007; Seow, 2013), further chemical reaction kinetics, and physical mechanisms like myosin head attachment (Nakajima et al., 1997) or ion diffusion, see the next sections 5.2 and 5.3. When trying to understand which process is effected by P_i kinetics, it is required to implement it in an even further enhanced model of muscle activation-contraction dynamics. The non-linear interactions of essential processes and mechanisms can only be understood in a cause-effect sense by mathematically formulating them in terms of at least algebraic but usually even differential equations and coupling these to state-of-the-art models formulated the same way. Eventually, this also applies to the very core of Hill-type models of muscle contraction: Hill's (phenomenological, steady-state) relation must then be replaced by a mechanistic model that represents both the basic force interactions and the thermodynamics on the cross-bridge level, formulated in terms of structure-based physical properties (cf. Günther et al., 2018).

5.2. Toward Resolving Characteristic Time Constants of Phosphate Dynamics

ODEs are used as standard method to model a characteristic evolution of an independent variable, which also depends on its own derivatives. Often, the time evolution of such an independent variable is studied in search for characteristic time constants of the modeled systems. The main Equation (3) of this contribution represents the time evolution of [ATP] or [P_i] with the characteristic time constants for condensation k_con and hydrolysis k_hyd. This ansatz together with the introduction of the state ${\tilde{\mu }}_{\text{ATP}}$ leads to the observed force decay $\frac{\text{d}}{\text{d}t}{F}_{MTU}$ (cf. Figures 1, 4). As described in the methods (section 3), the constants themselves describe a system characteristic, which represents a snapshot of even more detailed underlying processes.

Typical underlying processes for ATP hydrolysis k_hyd are the filling and draining of the ATP reservoirs including the actual type of flow during filling and draining: laminar or turbulent, the mixing of materials within the reservoir, and the residence time of the ATP on the myosin, i.e., the time ATP spends at the myosin head. Additionally, recent discoveries hint at two different types of ATP-myosin binding (Amrute-Nayak et al., 2014). The authors proposed that each myosin head has one site for ATP switching between two conformers (Tesi et al., 2017), which affects the binding duration and thus the characteristic time constant in the ODEs.

Admittedly, our model provides a rather phenomenological description of phosphate-induced force decay. In our activation model enhancement, the parameter ${\widehat{k}}_{\text{hyd}}$ determines the corresponding characteristic time of decay. It is remarkable that ${\widehat{k}}_{\text{hyd}}$ was the only free parameter in the enhancement part, that is, just its value was left to the parameter estimation process, whereas the other eight parameter values (cf. Table 1), were a priori fixed according to literature. The characteristic time 1/${\widehat{k}}_{\text{hyd}}$ of force decay is of the order 0.5 s. Force decay along with phosphate accumulation in the sarcoplasma is surely the physiological process that has a functional effect in natural contraction conditions. Yet, phosphate dilution can have the reciprocal impact of increasing isometric force (Phillips et al., 1993; Tesi et al., 2002). We suggest therefore to neutrally term the model enhancement “phosphate dynamics” rather than anything like “phosphate fatigue” or “phosphate-induced force decay.” The phosphate dynamics examined here is slow as compared to all other processes in the cross-bridge (Lymn-Taylor) cycle. Phosphate dynamics is thus rather a boundary condition for the Lymn-Taylor cycle than strongly and non-linearly interacting with cross-bridge dynamics. Therefore, the choice of our methodological path (see section 5.1) is a non-critical issue for our findings.

Other chemical factors that are known to have an effect on cross-bridge cycling, such as magnesium contributions, calcium precipitation, or glycolysis, have so far not been incorporated into our enhanced model of activation dynamics. Accordingly, the narrowing down to adding a single ODE for [ATP] dynamics seems to be nothing else than mathematically formulating the Lymn-Taylor cycle (Lymn and Taylor, 1971). Figure 5 gives a schematic overview of this cycle in accordance with our current interpretations (cf. processes P6 to P8 in Rockenfeller and Günther, 2016, Appendix A). Additionally, we included the illustration of a recent mechanistic idea of how the power stroke may be driven by Coulomb repulsion of ADP²⁻ and ${\text{P}}_i^{2-}$ (Rosenfeld, 2012; Günther et al., 2018). Our scheme is further consistent with crystallography measurements of geometric configurations (conformational states) of the myosin S1 part (Geeves and Holmes, 1999, pp. 700ff):

0: Without calcium available, tropomyosin blocks the access of myosin to new actin sites. Hence, myosin heads might be bound in a rigor formation, cf. state 6, or detached with either ATP (state 7) or its hydrolysis products present in the catalytic domain of the S1 region. In the former case, the switch-2 region of myosin is in an OPEN state, which facilitates phosphate release; in the latter case, the switch-2 region is CLOSED, i.e., forms some kind of phosphate-“pocket,” which enables ATP hydrolysis (Geeves and Holmes, 1999).

1: Soon after a calcium ion has bound to a troponin C terminal, several active sites on actin are exposed, where the exact number as well as the underlying mechanism are yet controversial, (Reconditi, 2006; Deasi et al., 2015). The gray continuation to the left of the actin filament accounts for the post-power-stroke translation after one full cycle (1-…-7), i.e., re-entering step 1 from step 7.

2: The catalytic domain binds to an available active site and changes its state from CLOSED to OPEN.

3: According to Elliott and Worthington (1994), Lampinen and Noponen (2005), and Rosenfeld (2012), the basic force that drives the power stroke is of electro-static (Coulomb) character: the OPEN state then allows the negatively charged phosphate to push away from the likewise negative charged ADP, and the Coulomb force is levered by the light chain domain to cause elastical deformations (bending) in the light chain domain itself as well as in myosin and actin filaments. All of these deformations summing up to ~4 nm pre-strain (Piazzesi and Lombardi, 1995; Piazzesi et al., 2002) in the isometric condition.

4: Eventually, as P_i and ADP move further away from each other, the actual power stroke (including release of the elastically stored energy) causes the S1 region to rotate relative to the S2 region and consequently moves the myosin rod by an equivalent of another about 7 nm (Günther et al., 2018, Figure 2) at its tip (attachment point to actin), with the whole power stroke length being altogether about 11 nm (Piazzesi and Lombardi, 1995; Piazzesi et al., 2002).

5: After having transformed the chemical free energy to mechanical work, phosphate is finally released to the sarcolemma (Takagi et al., 2004; Muretta et al., 2015). This process might be slowed down or hindered by higher [P_i] in the surrounding (Tesi et al., 2002).

6: In a consecutive step, ADP is likewise released into the sarcolemma (Suzuki et al., 1998), possibly with the aid of magnesium (Geeves and Holmes, 1999), leaving the actomyosin complex in a rigor formation.

7: In the presence of ATP, the myosin head detaches, switching again from OPEN to CLOSED state, allowing for another hydrolysis and the beginning of a new cycle. Just as step 5, this process can also be interfered with by a higher [P_i] (Kerrick and Xu, 2004). Note that ATP hydrolysis is commonly assumed to take place before re-attachment, but might also take place afterwards (Tonomura et al., 1966; Adelstein, 1980).

In all, the presented model (Equation 3) subsumes several even more detailed processes by assuming k_hyd and k_con. It remains open to integrate these into the equation, additionally. As was shown in this contribution, prerequisites are data of very good quality and a first guess of the system dynamics. Until then, the presented approach integrates phenomenologically-based molecular kinetics into macroscopic muscle models, enhancing them tremendously.

5.3. On the Effects of Phosphate Accumulation and Myofibrillar Calcium Sensitivity

The process of fatigue is commonly divided in three phases: An early decay of force, a plateau phase, and a late decay of force (Lännergren and Westerblad, 1991; Westerblad and Allen, 1991). The presented model considers the simulation of the early phase of fatigue. As the main cause for early force decay, single fiber experiments revealed (1) a decreased ability of the actomyosin crossbridges to generate force and (2) reduced myofibrillar [Ca²⁺] sensitivity (Allen et al., 2008). The first component is associated with elevated [P_i] due to the breakdown of PCr. When [P_i] had been elevated from below 1 to ~14 mM, a force decrease of about 50% was observed in [Ca²⁺]-activated skinned fiber experiments (Millar and Homsher, 1990, Figure 2). However, early skinned fiber experiments were conducted at non-physiological temperatures (10 to 15°C). The effect of elevated [P_i] on force might decrease with increasing temperature to physiological relevant conditions (Debold et al., 2004; Ranatunga, 2010). Nocella et al. (2017) examined the effect of [P_i] on the cross-bridge kinetics at physiological temperatures (33°C) and observed a fiber force decrease of 30% after 25 tetanic stimuli (on:off cycle = 0.4:1.5 s). Focusing on the time course of force decay (Nocella et al., 2017, Figure 2) after two tetanic stimuli (0.8 ms stimulation time, i.e., comparable to our presented experiments) they observed a force decrease of about 5%, which is similar to our observations (cf. Figure 1). Furthermore, they showed that the decrease in tetanic force mainly results from depressing the individual cross-bridge force and accelerated cross-bridge kinetics. However, force reduction during the early phase of fatigue was also associated with reduced myofibrillar [Ca²⁺] sensitivity (Debold, 2016). It was shown that the force decayed by 10% while the myoplasmic [Ca²⁺] increased (Westerblad and Allen, 1991). This increase is interpreted as the result of reduced myoplasmic [Ca²⁺] buffering (Westerblad and Allen, 1993). Alterations of the myofibrillar [Ca²⁺] sensitivity might occur due to increase of elevated metabolites as [H⁺] and [P_i]. Nelson and Fitts (2014) observed at a pH of 6.2 in skinned muscle fibers (at 30°C) a decreased myofibrillar [Ca²⁺] sensitivity. The myofibrillar [Ca²⁺] sensitivity decreased even more when they added 30 mM [P_i]. It was concluded that both low pH and elevated [P_i] have a substantial effect on myofibrillar [Ca²⁺] sensitivity. However, experiments with intact single fibers show a minor effect of acidose on tetanic force decrease (<10%) (Westerblad et al., 1997). The role of acidosis in acute fatigue remains controversial and a major unresolved issue is whether the force-reducing effects of elevated [P_i] in fatigue are amplified by the concomitant acidosis (Cheng et al., 2018). In our experiments, we can almost exclude the possibility that the pH decreased to values of 6.2 within a time period of 0.8 s (Stutzig et al., 2017). Thus, it seems that the main cause of fatigue that we observed and simulated is based on elevated [P_i], which influences both actomyosin cross-bridge force generation and myofibrillar [Ca²⁺] sensitivity.

5.4. Comparison of Muscle Parameters

Muscle parameters of the present study were determined by fitting the muscle model to a series of isometric contractions (n = 6…7). Typically a much higher number of experiments (n = 20…40, Scott et al., 1996; Curtin et al., 1998; Wagner et al., 2005; Siebert et al., 2008) must be used to determine these muscle model parameters, including for example isometric, isotonic, isokinetic, and quick-release experiments. Thus, from an experimental point of view, the herein applied model-based parameter estimation (fitting method) (Wagner et al., 2005; Siebert et al., 2007; Rockenfeller and Günther, 2016) is more efficient compared to the classic procedure. Rabbit GAS and PLA muscle parameters have been previously determined with classic methods (Siebert et al., 2015). In general, their results are in good agreement with muscle model parameters determined in the present study. Maximum shortening velocity was overestimated (GAS: factor 1.7; PLA: factor 3.1) compared to classic methods, but consistent with recent findings on non-steady-state contractions (Piazzesi et al., 2002; Rockenfeller and Günther, 2016; Günther et al., 2018). Furthermore, maximum power deviates by 25% compared to classic methods. This can be explained by the limited parameter range available in the fitting method. Shortening velocity of the contractile component reaches maximally 0.5 v_max in isometric contractions used for parameter fitting. Thus, uncertainty in parameter estimation is higher for v_max and P_max compared to classic methods, in which v_max can be approximated by contractions against low loads or unloaded contractions (Edman, 1979). In contrast, parameters of the force-length relation are almost similar between classic and fitting method. There were differences in optimum muscle length (GAS: 11%; PLA: 20%) and widths of the ascending limb (GAS: 2%, PLA: 17%). For the SEE characteristic, our fit yielded a rather broad non-linear toe region (ΔU_{SEE, l} ≈ 8 vs. 4.9% for GAS and 3.6% for PLA with classic methods), with a transition to the linear region at around F_max and a less stiff behavior in the linear region (K_{SEE, l} ≈ 22.4 kN/m Günther et al., 2007, Equation 4 vs. k ≈ 30.3 kN/m Siebert et al., 2015, Equation 5 for GAS and K_{SEE, l} ≈ 15.6 kN/m vs. k ≈ 21.9 kN/m for PLA). Determination of PEE characteristics is limited to the passive force range (GAS: 20% F_max; PLA: 9% F_max) covered by isometric experiments at different muscle lengths (see Figure 1). Differences in PEE stiffness at longest muscle length is 10% for GAS and 25% for PLA. Deviations in fitted compared to experimental (classic method) PEE stiffness of PLA occurred due to lower passive force range available for parameter fitting (compared to GAS) as well as comparably low PEE forces (<5 N) and thus small impact on model simulations. Thus, there are only small differences in PEE forces (equaling passive force in the inactive muscle) between experiment and model simulation, see Figure 1 at t < 0.

5.5. Length-Dependence of Fatigue

Muscular fatigue plays an important role in the assessment of work-place ergonomics in order to accurately predict demands on workers with respect to the muscular forces required for their work tasks. To this end, endurance time is a measure used to characterize muscular loading situations. It was introduced to quantify the time a subject can hold a specific load by muscular contraction (Rohmert, 1960) and has been used to study many different postures and muscles (e.g., Frey-Law and Avin, 2010). In general, measurements reveal that endurance time is shorter for higher muscular forces. Furthermore, experiments show that below a certain load, endurance time becomes very long indicating a muscular load where the normal ATP resynthesis rate is sufficient to compensate the static energy requirement for the muscles. This lower threshold may be somewhere between 2 and 20% of the maximum voluntary contraction (van Dieën and oude Vrielink, 1994; El ahrache et al., 2006). The model presented here also shows this behavior (Figure 6), although we prefer the term exhaustion time when talking about the inability to maintain a certain, length-dependent force. Herein, exhaustion time was defined to be the first time instant where force had decayed more than 5% of its initial value at t = 0 s. For stimulation values around 0.1, exhaustion time may rise well above 20 s.

FIGURE 6

Figure 6. Simulated time courses for CE length (A) and MTU force (B) for a GAS muscle during isometric contractions, utilizing fixed parameters from Table 2. CE length are displayed relative to ℓ_{CE, opt} ≈ 0.0198m and forces relative to F_max ≈ 133N. Muscle length varies about −10, …, +15mm, in steps of 5 mm, around the reference length ℓ_MTU = 122mm. Colors and line-styles in both sub-figures are coherent, indicating MTU lengths [as specified in (A)] and stimulation levels [as specified in (B)]. Two observations should be highlighted. First, the force level in (B) at which (theoretically) infinite exhaustion time occurs settles at around 0.16, corresponding to 16%F_max, which is in perfect agreement with the 15%F_max from in vivo experimental data (Rohmert, 1960). Second, the force at the longest muscle length [blue line in (B)] settles at significantly higher values, which is due to passive (PEE) forces that were herein not modeled to show any exhaustion.

In addition, the model predicts a length dependence of the exhaustion time (Figure 7). One core assumption of the model is that the hydrolysis rate decreases with increasing muscle fiber length (see Equation 4). This assumption was derived from the experimental observation that the force decay due to exhaustion scales with muscle length, see experimental data in Figure 4. This characteristic is, thus, immediately reproduced by the model, as shown in Figure 4 and Equation (8). More precisely, our model strongly suggests that the mechanism(s) that govern the exhaustion process are the same as those governing the length dependency of ATP hydrolysis.

FIGURE 7

Figure 7. (A) Simulated force-length-exhaustion time diagram for variations in MTU lengths and stimulation levels as given in Figure 6. (B–D) Contour plots through the three different planes in (A). In (B), the force-length diagram for various exhaustion times is shown. The active and passive force-length characteristic of the CE is clearly visible. In (C), the force-exhaustion time courses for various CE length are displayed. The longer the CE, the longer the exhaustion time. The typical exponential (or hyperbolic) characteristic (Rohmert, 1960; Frey-Law and Avin, 2010) is visible. (D) Shows CE length-exhaustion time curves for various force levels. The higher the force, the shorter the exhaustion time at the same length. For longer muscles, the passive force determines the exhaustion time.

This is an issue also discussed in the literature with respect to endurance time. Interestingly, the findings are ambiguous. Some experiments also found that fatiguability is reduced for longer muscle lengths, i.e., endurance time increases with muscle length (in accordance with our model, cf. Figure 7D) (Sacco et al., 1985; de Haan et al., 1986; Matthijsse et al., 1987; McKenzie and Gandevia, 1987; Arendt-Nielsen et al., 1992; Kawakami et al., 2000; Rassier, 2000). However, a few experiments also indicate that endurance time is reduced for increasing muscle length (Fitch and McComas, 1985; Ng et al., 1985; Willems and Stauber, 2002; Kooistra et al., 2005; Lee et al., 2007). This would be in contrast to our model prediction. In many of these studies, however, only few lengths were investigated omitting a direct conclusion on the length dependence. Finally, the data shown in Petrofsky and Phillips (1980) for endurance time in relation to elbow angle indicates that there is a maximum endurance time for an elbow angle of 120°, which is also the elbow angle for which the subjects strength is maximal. This indicates, that the endurance time may be related to the force-length relation of the muscle, an effect also visible in our model.

There is one important difference though between in-vivo endurance time in human workers and the exhaustion time predicted by our model: While muscular fatigue in humans may (partially) be compensated for by increasing motor-unit recruitment, i.e., increased muscle stimulation, we here simulate exhaustion under constant muscle stimulation input (u = 1). This additional recruitment is crucial for the notion of endurance time (Petrofsky, 1978) and can be seen in humans using muscle surface electromyograms (EMG) (Gamet and Maton, 1989; Maton and Gamet, 1989). It is actually used as an indicator of fatigue in ergonomics (Jørgensen, 1997; Nussbaum et al., 2001; Garg et al., 2002). While this aspect has to be considered in ergonomics, our model approach, for the first time, allows a systematic model-based analysis of the length-dependence of muscular exhaustion.

6. Summary

In this work, we developed a model able to explain the time course of force decay, which occurs as a consequence of ongoing neural stimulation. As opposed to the widespread but rather diffuse term fatigue, we here prefer to term the muscle's incapability of maintaining a certain force level as exhaustion in case the dominating mechanism behind fatigue is a shift in the equilibrium of the ATP hydrolysis-condensation kinetics (phosphate dynamics). Accordingly, we assume as the key feature of exhaustion the deteriorating chemical potential, i.e. — the change in internal energy per change in particle number — of ATP during hydrolysis. More precisely, decreasing [ATP] as well as increasing either [ADP] or [P_i] yields a lower chemical potential being reflected in exhaustion. We incorporated phosphate dynamics into an established Hill-type muscle model representing excitation-contraction dynamics. This new (merged and enhanced) model was validated by parameter estimation with using experimental data of isometric contractions gathered from two types of rabbit calf muscles. With parameter values obtained from optimally fitting direct dynamic model simulations to experimental contraction data, the model can reproduce experimental findings strikingly better than the initial model. Moreover, the parameter values well agree with prior estimates from literature and eventually allow for predicting measurements from experiments with longer stimulation duration. We argue that the presented methodology of model enhancement can and ought to be applied to further physiological mechanisms, e.g., the Lymn-Taylor cycle or the impact of changes in myofibrillar calcium concentration. Lastly, consequences of the length-dependencies within the model have been investigated and linked to known findings from ergonomics.

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Ethics Statement

Experiments were approved according to section 8 German animal protection law (Tierschutzgesetz, BGBlI 1972, 1277) by the competent authority for animal welfare in Thuringia, Germany (Landesamt fur Verbraucherschutz, Abteilung Gesundheitlicher und technischer Verbraucherschutz; Permit Number: 02-022/11 and 02-027/14).

Author Contributions

RR developed the model idea as well as the first draft of the manuscript and coordinated the partners. MG contributed the physiological insights, particularly during the parameter optimization and within section 5.1. DH investigated the aspect of length-dependency within the model as well as the inter-linkage with ergonomics, particularly in section 5.5. NS discussed the biological ramifications (cf. section 5.3) thus assuring to keep the model as straightforward as possible. TS, KL, and MB provided the data. TS further assessed the optimized parameter values (see section 5.4). SS provided the valuable discussions, final revision, and co-authored section 5.2. KL conducted the experiments and supervised the methodological aspects. MB conducted the several discussions as well as a final revision. TG helped in analyzing the mathematical methods (see the Appendices).

Funding

MG was supported by the Cluster of Excellence for Simulation Technology (SimTech) (DFG: EXC310), Deutsche Gesetzliche Unfallversicherung (DGUV) project Wirbelsäulenmodell passive Strukturen, and Deutsche Forschungsgemeinschaft (DFG: SCHM2392/5-2), all granted to SS. DH was supported by the Ministerium für Wissenschaft, Forschung und Kunst Baden-Württemberg (MWK, Az: 33-7533.-30-20/7/2). TS was supported by DFG grants SI841/3-2 and SI841/7-1,2. MB was supported by DFG grants BO3091/4-1 and BO3091/4-2. SS was funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2075 - 390740016. The publication costs were funded by the Open Access Fund of the University of Koblenz-Landau.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Abbiss, C. R., and Laursen, P. B. (2005). Models to explain fatigue during prolonged endurance cycling. Sports Med. 35, 865–898. doi: 10.2165/00007256-200535100-00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Adelstein, R. S. (1980). Regulation and kinetics of the actin-myosin-ATP interaction. Annu. Rev. Biochem. 49, 921–956. doi: 10.1146/annurev.bi.49.070180.004421

PubMed Abstract | CrossRef Full Text | Google Scholar

Aljure, E. F., and Borrero, L. M. (1968). The influence of muscle length on the development of fatigue in toad sartorius. J. Physiol. 199, 241–252. doi: 10.1113/jphysiol.1968.sp008651

CrossRef Full Text | Google Scholar

Allen, D. G., Lamb, G. D., and Westerblad, H. (2008). Skeletal muscle fatigue: cellular mechanisms. Physiol. Rev. 88, 287–332. doi: 10.1152/physrev.00015.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, D. G., and Orchard, C. H. (1987). Myocardial contractile function during ischemia and hypoxia. Circ. Res. 60, 153–168. doi: 10.1161/01.RES.60.2.153

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, D. G., and Westerblad, H. (2001). Role of phosphate and calcium stores in muscle fatigue. J. Physiol. 536, 657–665. doi: 10.1111/j.1469-7793.2001.t01-1-00657.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Amrute-Nayak, M., Lambeck, K.-A., Radocaj, A., Huhnt, H. E., Scholz, T., Hahn, N., et al. (2014). Atp turnover by individual myosin molecules hints at two conformers of the myosin active site. Proc. Natl. Acad. Sci. U.S.A. 111, 2536–2541. doi: 10.1073/pnas.1316390111

PubMed Abstract | CrossRef Full Text | Google Scholar

Arendt-Nielsen, L., Gantchev, N., and Sinkær, T. (1992). The influence of muscle length on muscle fibre conduction velocity and development of muscle fatigue. Electroencephalogr. Clin. Neurophysiol. 85, 166–172. doi: 10.1016/0168-5597(92)90128-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Azizi, E., Brainerd, E. L., and Roberts, T. J. (2008). Variable gearing in pennate muscles. Proc. Natl. Acad. Sci. U.S.A. 105, 1745–1750. doi: 10.1073/pnas.0709212105

PubMed Abstract | CrossRef Full Text | Google Scholar

Bagshaw, C. R., and Trentham, D. R. (1973). The reversibility of adenosine triphosphate cleavage by myosin. Biochem. J. 133, 323–328. doi: 10.1042/bj1330323

PubMed Abstract | CrossRef Full Text | Google Scholar

Barclay, C. J. (2015). Energetics of contraction. Compr. Physiol. 5, 961–995. doi: 10.1002/cphy.c140038

PubMed Abstract | CrossRef Full Text | Google Scholar

Baylor, S. M., and Hollingworth, S. (1998). Model of sarcomeric Ca²⁺ movements, including ATP Ca²⁺ binding and diffusion, during activation of frog skeletal muscle. J. Gen. Physiol. 112, 297–316. doi: 10.1085/jgp.112.3.297

PubMed Abstract | CrossRef Full Text | Google Scholar

Bechtel, P. J., and Best, P. M. (1985). “Integration of ATP production and consumption,” in Reciprocal Meat Conference Proceedings (Baton Rouge, LA), Vol. 38, 26–31.

Google Scholar

Bergström, M., and Hultman, E. (1985). Energy cost and fatigue during intermittent electrical stimulation of human skeletal muscle. J. Appl. Physiol. 65, 1500–1505. doi: 10.1152/jappl.1988.65.4.1500

PubMed Abstract | CrossRef Full Text | Google Scholar

Bickham, D. C. (2003). Extracellular K⁺ accumulation: a physiological framework for fatigue during intense exercise. J. Physiol. 554:593. doi: 10.1113/jphysiol.2003.059139

PubMed Abstract | CrossRef Full Text | Google Scholar

Bickham, D. C., West, T. G., Webb, M. R., Woledge, R. C., Curtin, N. A., and Ferenczi, M. A. (2011). Milisecond-scale biochemical response to change in strain. Biophys. J. 101, 2445–2454. doi: 10.1016/j.bpj.2011.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bigland, B., and Lippold, O. C. J. (1954). Motor unit activity in the voluntary contraction of human muscle. J. Physiol. 125, 322–335. doi: 10.1113/jphysiol.1954.sp005161

PubMed Abstract | CrossRef Full Text | Google Scholar

Blümel, M., Hooper, S. L., Guschlbauer, C., White, W. W., and Büschges, A. (2012). Determining all parameters necessary to build Hill-type muscle models from experiments on single muscles. Biol. Cybern. 106, 543–558. doi: 10.1007/s00422-012-0531-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M. (2009). A finite element approach for the modelling of skeletal muscle fatigue. GAMM Mitteil. 32, 205–220. doi: 10.1002/gamm.200910017

CrossRef Full Text | Google Scholar

Böl, M., Leichsenring, K., Ernst, M., Wick, C., Blickhan, R., and Siebert, T. (2015). Novel microstructural findings in M. plantaris and their impact during active and passive loading at the macro level. J. Mech. Behav. Biomed. Mater. 51, 25–39. doi: 10.1016/j.jmbbm.2015.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M., Leichsenring, K., Weichert, C., Sturmat, M., Schenk, P., Blickhan, R., et al. (2013). Three-dimensional surface geometries of the rabbit soleus muscle during contraction: input for biomechanical modelling and its validation. Biomech. Model. Mechanobiol. 12, 1205–1220. doi: 10.1007/s10237-013-0476-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, I. E., and Loeb, G. E. (2000). Measured and modeled properties of mammalian skeletal muscle: IV. Dynamics of activation and deactivation. J. Muscle Res. Cell Motil. 21, 33–47. doi: 10.1023/A:1005687416896

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, N., Bubeck, D., Haeufle, D. F., Weickenmeier, J., Kuhl, E., Alt, W., et al. (2017). Weekly time course of neuro-muscular adaptation to intensive strength training. Front. Physiol. 8:329. doi: 10.3389/fphys.2017.00329

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, Z., Bai, E., and Shields, R. K. (2009). Fatigue and non-fatigue mathematical muscle models during functional electrical stimulation of paralyzed muscle. Biomed. Signal Process. Control 5, 87–93. doi: 10.3182/20090812-3-DK-2006.0007

PubMed Abstract | CrossRef Full Text | Google Scholar

Callahan, D. M., Umberger, B. R., and Kent, J. A. (2016). Mechanisms of in vivo muscle fatigue in humans: investigating age-related fatigue resistance with a computational model. J. Physiol. 594, 3407–3421. doi: 10.1113/JP271400

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardozo, A. C., Gonçalves, M., and Dolan, P. (2011). Back extensor muscle fatigue at submaximal workloads assessed using frequency banding of the electromyographic signal. Clin. Biomech. 26, 971–976. doi: 10.1016/j.clinbiomech.2011.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlson, F. D., and Wilkie, D. R. (1974). Muscle Physiology. Englewood Cliffs, NJ: Prentice-Hall Inc.

Google Scholar

Chen, S. K., Simonian, P. T., Wickiewicz, T. L., Otis, J. C., and Warren, R. F. (1999). Radiographic evaluation of glenohumeral kinematics: a muscle fatigue model. J. Elbow Shoulder Surg. 8, 49–52. doi: 10.1016/S1058-2746(99)90055-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, A. J., Place, N., and Westerblad, H. (2018). Molecular basis for exercise-induced fatigue: The importance of strictly controlled cellular Ca²⁺ handling. Cold Spring Harbor Perspect. Med. 8:a029710. doi: 10.1101/cshperspect.a029710

PubMed Abstract | CrossRef Full Text | Google Scholar

Cooke, R. (2007). Modulation of the actomyosin interaction during fatigue of skeletal muscle. Muscle Nerve 36, 756–777. doi: 10.1002/mus.20891

PubMed Abstract | CrossRef Full Text | Google Scholar

Curtin, N. A., Gardner-Medwin, A. R., and Woledge, R. C. (1998). Predictions of the time course of force and power output by dogfish white muscle fibres during brief tetani. J. Exp. Biol. 201, 103–114.

PubMed Abstract | Google Scholar

de Haan, A., de Jong, J., van Doorn, J. E., Huijing, P. A., Woittiez, R. D., and Westra, H. G. (1986). Muscle economy of isometric contractions as a function of stimulation time and relative muscle length. Pflügers Archiv. Eur. J. Physiol. 407, 445–450. doi: 10.1007/BF00652632

PubMed Abstract | CrossRef Full Text | Google Scholar

De La Cruz, E. M., and Ostap, E. M. (2009). Kinetic and equilibrium analysis of the myosin ATPase. Methods Enzymol. 455, 157–192. doi: 10.1016/S0076-6879(08)04206-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Deasi, R., Geeves, M. A., and Kad, N. M. (2015). Using fluorescent myosin to directly visualize cooperative activation of thin filaments. J. Biol. Chem. 290, 1915–1925. doi: 10.1074/jbc.M114.609743

CrossRef Full Text | Google Scholar

Debold, E. P. (2016). Decreased myofilament calcium sensitivity plays a significant role in muscle fatigue. Exerc. Sport Sci. Rev. 44, 144–149. doi: 10.1249/JES.0000000000000089

PubMed Abstract | CrossRef Full Text | Google Scholar

Debold, E. P., Dave, H., and Fitts, R. H. (2004). Fiber type and temperature dependence of inorganic phosphate: implications for fatigue. Am. J. Cell Physiol. 287, C673–C681. doi: 10.1152/ajpcell.00044.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Debold, E. P., Fitts, R. H., Sundberg, C. W., and Nosek, T. M. (2016). Muscle fatigue from the perspective of a single crossbridge. Med. Sci. Sports Exerc. 48, 2270–2280. doi: 10.1249/MSS.0000000000001047

PubMed Abstract | CrossRef Full Text | Google Scholar

Deeb, J. M., Drury, C. G., and Pendergast, D. R. (1992). An exponential model of isometric muscular fatigue as a function of age and muscle groups. Ergonomics 35, 899–918. doi: 10.1080/00140139208967370

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., Wexler, A. S., and Binder-MacLeod, S. A. (2000). A predictive model of fatigue in human skeletal muscles. J. Appl. Physiol. 89, 1322–1332. doi: 10.1152/jappl.2000.89.4.1322

PubMed Abstract | CrossRef Full Text | Google Scholar

Dorgan, S. J., and O'Malley, M. J. (1998). A mathematical model for skeletal muscle activated by n-let pulse trains. IEEE Trans. Rehabil. Eng. 6, 286–299. doi: 10.1109/86.712226

PubMed Abstract | CrossRef Full Text | Google Scholar

Doud, J. R., and Walsh, J. M. (1995). Muscle fatigue and length interaction: effect on the EMG frequency component. Electromyogr. Clin. Neurophysiol. 35, 331–339.

PubMed Abstract | Google Scholar

Dragomir, C. T. (1970). On the nature of forces acting between myofilaments in resting state and under contraction. J. Theor. Biol. 27, 343–356. doi: 10.1016/S0022-5193(70)80001-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Drazan, J. F., Hullfish, T. J., and Baxter, J. R. (2019). An automatic fascicle tracking algorithm quantifying gastrocnemius architecture during maximal effort contractions. PeerJ 7:e7121. doi: 10.7717/peerj.7120

PubMed Abstract | CrossRef Full Text | Google Scholar

Edman, K. A. P. (1979). The velocity of unloaded shortening and its relation to sarcomere length and isometric force in vertebrate muscle fibres. J. Physiol. 291, 143–159. doi: 10.1113/jphysiol.1979.sp012804

PubMed Abstract | CrossRef Full Text | Google Scholar

El ahrache, K., Imbeau, D., and Farbos, B. (2006). Percentile values for determining maximum endurance times for static muscular work. Int. J. Ind. Ergon. 36, 99–108. doi: 10.1016/j.ergon.2005.08.003

CrossRef Full Text | Google Scholar

Elliott, G., and Worthington, C. (1994). How muscle may contract. Biochim. Biophys. Acta 1200, 109–116. doi: 10.1016/0304-4165(94)90124-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Enoka, R. M., and Stuart, D. G. (1992). Neurobiology of muscle fatigue. J. Appl. Physiol. 72, 1631–1648. doi: 10.1152/jappl.1992.72.5.1631

PubMed Abstract | CrossRef Full Text | Google Scholar

Eriksson, A., Holmberg, H. C., and Westerblad, H. (2016). A numerical model for fatigue effects in whole-body human exercise. Math. Comput. Modell. Dyn. Syst. 22, 21–38. doi: 10.1080/13873954.2015.1083592

CrossRef Full Text | Google Scholar

Fenwick, A. J., Leighton, S. R., and Tanner, B. C. W. (2016). Myosin MgADP release rate decreases as sarcomere length increases in skinned rat soleus muscle fibers. Biophys. J. 111, 2011–2023. doi: 10.1016/j.bpj.2016.09.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitch, S., and McComas, A. (1985). Influence of human muscle length on fatigue. J. Physiol. 362, 205–213. doi: 10.1113/jphysiol.1985.sp015671

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitts, R. H. (1994). Cellular mechanisms of muscle fatigue. Physiol. Rev. 74, 49–94. doi: 10.1152/physrev.1994.74.1.49

PubMed Abstract | CrossRef Full Text | Google Scholar

Freund, J., and Takala, E. (2001). A dynamic model of the forearm including fatigue. J. Biomech. 34, 597–605. doi: 10.1016/S0021-9290(01)00009-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Frey-Law, L. A., and Avin, K. G. (2010). Endurance time is joint-specific: a modelling and meta-analysis investigation. Ergonomics 53, 109–129. doi: 10.1080/00140130903389068

PubMed Abstract | CrossRef Full Text | Google Scholar

Frey-Law, L. A., Looft, J. M., and Heitsman, J. (2012). A three-compartment muscle fatigue model accurately predicts joint-specific maximum endurance times for sustained isometric tasks. J. Biomech. 45, 1803–1808. doi: 10.1016/j.jbiomech.2012.04.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Gadian, D. G., Radda, G. K., Brown, T. R., Chance, E. M., Dawson, M. J., and Wilkie, D. R. (1981). The activity of creatine kinase in frog skeletal muscle studied by saturation-transfer nuclear magnetic resonance. Biochem. J. 194, 215–228. doi: 10.1042/bj1940215

PubMed Abstract | CrossRef Full Text | Google Scholar

Gamet, D., and Maton, B. (1989). The fatigability of two agonistic muscles in human isometric voluntary submaximal contraction: an EMG study. I. Assessment of muscular fatigue by means of surface EMG. Eur. J. Appl. Physiol. Occupat. Physiol. 58, 361–368. doi: 10.1007/BF00643510

PubMed Abstract | CrossRef Full Text | Google Scholar

Gandevia, S. (2001). Spinal and supraspinal factors in human muscle fatigue. Physiol. Rev. 81, 1725–1789. doi: 10.1152/physrev.2001.81.4.1725

PubMed Abstract | CrossRef Full Text | Google Scholar

Gandevia, S. C., Enoka, R. M., McComas, A. J., Stuart, D. G., Thomas, C. K., Pierce, P. A., et al. (1995). Fatigue: Neural and Muscular Mechanisms. New York, NY: Springer US.

Google Scholar

Garg, A., Hegmann, K. T., Schwoerer, B. J., and Kapellusch, J. M. (2002). The effect of maximum voluntary contraction on endurance times for the shoulder girdle. Int. J. Ind. Ergon. 30, 103–113. doi: 10.1016/S0169-8141(02)00078-1

CrossRef Full Text | Google Scholar

Geeves, M. A., and Holmes, K. C. (1999). Structural mechanism of muscle contraction. Annu. Rev. Biochem. 68, 687–728. doi: 10.1146/annurev.biochem.68.1.687

PubMed Abstract | CrossRef Full Text | Google Scholar

Giat, Y., Mizrahi, J., and Levy, M. (1993). A musculotendon model of the fatigue profiles of paralyzed quadriceps muscle under FES. IEEE Trans. Biomed. Eng. 40, 664–674. doi: 10.1109/10.237696

PubMed Abstract | CrossRef Full Text | Google Scholar

Glass, L. D. (2017). The role of calcium in the mechanism of muscle fatigue in intact, single muscle fibres of mice at physiological temperature (Ph.D. thesis), University of Calgary, Calgary, AB, Canada.

Google Scholar

Grigorenko, B. L., Rogov, A. V., Topol, I. A., Burt, S. K., and Martinez, H. M. (2007). Mechanism of the myosin catalyzed hydrolysis of ATP as rationalized by molecular modeling. Proc. Natl. Acad. Sci. U.S.A. 104, 7057–7061. doi: 10.1073/pnas.0701727104

PubMed Abstract | CrossRef Full Text | Google Scholar

Günther, M., Haeufle, D. F. B., and Schmitt, S. (2018). The basic mechanical structure of the skeletal muscle machinery: one model for linking microscopic and macroscopic scales. J. Theor. Biol. 456, 137–167. doi: 10.1016/j.jtbi.2018.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Günther, M., and Schmitt, S. (2010). A macroscopic ansatz to deduce the Hill relation. J. Theor. Biol. 263, 407–418. doi: 10.1016/j.jtbi.2009.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Günther, M., Schmitt, S., and Wank, V. (2007). High-frequency oscillations as a consequence of neglected serial damping in Hill-type muscle models. Biol. Cybern. 97, 63–79. doi: 10.1007/s00422-007-0160-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Guynn, R. W., and Veech, R. L. (1973). The equilibrium constants of the adenosine triphosphate hydrolysis and the adenosine triphosphate-citrate lyase reactions. J. Biol. Chem. 248, 6966–6972.

PubMed Abstract | Google Scholar

Haeufle, D. F. B., Günther, M., Bayer, A., and Schmitt, S. (2014a). Hill-type muscle model with serial damping and eccentric force-velocity relation. J. Biomech. 47, 1531–1536. doi: 10.1016/j.jbiomech.2014.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Haeufle, D. F. B., Günther, M., Wunner, G., and Schmitt, S. (2014b). Quantifying control effort of biological and technical movements: an information-entropy-based approach. Phys. Rev. E 89:012716. doi: 10.1103/PhysRevE.89.012716

PubMed Abstract | CrossRef Full Text | Google Scholar

Hancock, C. R., Janssen, E., and Terjung, R. L. (2005). Skeletal muscle contractile performance and ADP accumulation in adenylate kinase-deficient mice. Am. J. Cell Physiol. 288, C1287–C1297. doi: 10.1152/ajpcell.00567.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatze, H. (1977). A myocybernetic control model of skeletal muscle. Biol. Cybern. 25, 103–119. doi: 10.1007/BF00337268

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatze, H. (1978). A general myocybernetic control model of skeletal muscle. Biol. Cybern. 28, 143–157. doi: 10.1007/BF00337136

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatze, H. (1981). Myocybernetic Control Models of Skeletal Muscle. Pretoria: University of South Africa.

Google Scholar

Hawkins, D., and Hull, M. L. (1993). Muscle forces as affected by fatigue: mathematical model and experimental verification. J. Biomech. 26, 1117–1128. doi: 10.1016/S0021-9290(05)80010-7

CrossRef Full Text | Google Scholar

Hibberd, M. G., Dantzig, J. A., Trentham, D. R., and Goldman, Y. E. (1985). Phosphate release and force generation in skeletal muscle fibers. Science 228, 1317–1319. doi: 10.1126/science.3159090

PubMed Abstract | CrossRef Full Text | Google Scholar

Hilber, K., Sun, Y. B., and Irving, M. (2001). Effects of sarcomere length and temperature on the rate of ATP utilisation by rabbit psoas muscle fibres. J. Physiol. 531, 771–780. doi: 10.1111/j.1469-7793.2001.0771h.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, A. V. (1910). The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. Proc. Physiol. Soc. 1, iv–vii.

Google Scholar

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. Proc. R. Soc. Lond. B 126, 136–195. doi: 10.1098/rspb.1938.0050

CrossRef Full Text | Google Scholar

Huxley, A. F. (1974). Muscular contraction. J. Physiol. 243, 1–43. doi: 10.1113/jphysiol.1974.sp010740

CrossRef Full Text | Google Scholar

James, A., and Green, S. (2012). A phenomenological model of muscle fatigue and the power-endurance relationship. J. Appl. Physiol. 113, 1643–1651. doi: 10.1152/japplphysiol.00800.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Jørgensen, K. (1997). Human trunk extensor muscles physiology and ergonomics. Acta Physiol. Scand. 637, 1–58.

PubMed Abstract | Google Scholar

Karatzaferi, C., de Haan, A., Ferguson, R. A., van Mechelen, W., and Sargeant, A. J. (2001). Phosphocreatine and ATP content in human single muscle fibres before and after maximum dynamic exercise content in human single muscle fibres before and after maximum dynamic exercise. Eur. J. Physiol. 442, 467–474. doi: 10.1007/s004240100600

PubMed Abstract | CrossRef Full Text | Google Scholar

Katz, B. (1939). The relation between force and speed in muscular contraction. J. Physiol. 96, 45–64. doi: 10.1113/jphysiol.1939.sp003756

PubMed Abstract | CrossRef Full Text | Google Scholar

Kawakami, Y., Amemiya, K., Kanehisa, H., Ikegawa, S., and Fukunaga, T. (2000). Fatigue responses of human triceps surae muscles during repetitive maximal isometric contractions. J. Appl. Physiol. 88, 1969–1975. doi: 10.1152/jappl.2000.88.6.1969

PubMed Abstract | CrossRef Full Text | Google Scholar

Kent-Braun, J. A. (1999). Central and peripheral contributions to muscle fatigue in humans during sustained maximal effort. Eur. J. Appl. Physiol. 80, 57–63. doi: 10.1007/s004210050558

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerrick, W. G. L., and Xu, Y. (2004). Inorganic phosphate affects the pCa-force relationship more than the pCa-ATPase by increasing the rate of dissociation of force generating cross-bridges in skinned fibers from both EDL and soleus muscles of the rat. J. Muscle Res. Cell Motil. 25, 107–117. doi: 10.1023/B:JURE.0000035841.04314.16

PubMed Abstract | CrossRef Full Text | Google Scholar

Kistemaker, D. A., van Soest, A. J., and Bobbert, M. F. (2005). Length-dependent [Ca²⁺] sensitivity adds stiffness to muscle. J. Biomech. 38, 1816–1821. doi: 10.1016/j.jbiomech.2004.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Komi, P. V. (1973). “Biomechanics III,” in 3rd International Seminar on Biomechanics, Volume 8, Chapter Measurement of the Force-Velocity Relationship in Human Muscle under Concentric and Eccentric Contractions (Basel: Karger), 224–229. doi: 10.1159/000393754

CrossRef Full Text | Google Scholar

Komi, P. V. (2000). Stretch-shortening cycle: a powerful model to study normal and fatigued muscle. J. Biomech. 33, 1197–1206. doi: 10.1016/S0021-9290(00)00064-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kooistra, R. D., de Ruiter, C. J., and de Haan, A. (2005). Muscle activation and blood flow do not explain the muscle length-dependent variation in quadriceps isometric endurance. J. Appl. Physiol. 98, 810–816. doi: 10.1152/japplphysiol.00712.2004

CrossRef Full Text | Google Scholar

Kowaltowski, A. J., Castilho, R. F., Grijalba, M. T., Bechara, E. J. H., and Vercesi, A. E. (1996). Effect of inorganic phosphate concentration on the nature of inner mitochondrial membrane alterations mediated by Ca²⁺ ions. J. Biol. Chem. 271, 2929–2934. doi: 10.1074/jbc.271.6.2929

PubMed Abstract | CrossRef Full Text | Google Scholar

Kübler, W., and Katz, A. M. (1977). Mechanism of early “pump” failure of the lschemic heart: possible role of adenosine triphosphate depletion and inorganic phosphate accumulation. Am. J. Cardiol. 40, 467–471. doi: 10.1016/0002-9149(77)90174-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lampinen, M. J., and Noponen, T. (2005). Electric dipole theory and thermodynamics of actomyosin molecular motor in muscle contraction. J. Theor. Biol. 236, 397–421. doi: 10.1016/j.jtbi.2005.03.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Lännergren, J., and Westerblad, H. (1991). Force decline due to fatigue and intracellular acidification in isolated fibres from mouse skeletal muscle. J. Physiol. 434, 307–322. doi: 10.1113/jphysiol.1991.sp018471

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawson, J. W. R., and Veech, R. L. (1979). Effects of pH and free Mg²⁺ on the K_eq of the creatine kinase reaction and other phosphate hydrolyses and phosphate transfer reactions. J. Biol. Chem. 254, 6528–6537.

PubMed Abstract | Google Scholar

Lee, S. C. K., Braim, A., Becker, C. N., Prosser, L. A., Tokay, A. M., and Binder-MacLeod, S. A. (2007). Diminished fatigue at reduced muscle length in human skeletal muscle. Muscle Nerve 36, 789–797. doi: 10.1002/mus.20873

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemaire, K. K., Baan, G. C., Jaspers, R. T., and van Soest, A. J. (2016). Comparison of the validity of Hill and Huxley muscle-tendon complex models using experimental data obtained from rat M. soleus in situ. J. Exp. Biol. 219, 977–987. doi: 10.1242/jeb.128280

CrossRef Full Text | Google Scholar

Linari, M., Caremani, M., and Lombardi, V. (2010). A kinetic model that explains the effect of inorganic phosphate on the mechanics and energetics of isometric contraction of fast skeletal muscle. Proc. R. Soc. B 277, 19–27. doi: 10.1098/rspb.2009.1498

PubMed Abstract | CrossRef Full Text | Google Scholar

Lindström, L., Kadefors, R., and Petersén, I. (1977). An electromyographic index for localized muscle fatigue. J. Appl. Physiol. 43, 750–754. doi: 10.1152/jappl.1977.43.4.750

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J. Z., Brown, R. W., and Yue, G. H. (2002). A dynamical model of muscle activation, fatigue, and recovery. Biophys. J. 82, 2344–2359. doi: 10.1016/S0006-3495(02)75580-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Lymn, R. W., and Taylor, E. W. (1971). Mechanism of adenosine triphosphate hydrolysis by actomyosin. Biochemistry 10, 4617–4624. doi: 10.1021/bi00801a004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, L., Chablat, D., Bennis, F., and Zhang, W. (2009). A new simple dynamic muscle fatigue model and its validation. Int. J. Ind. Ergon. 39, 211–220. doi: 10.1016/j.ergon.2008.04.004

CrossRef Full Text | Google Scholar

Ma, L., Chablat, D., Bennis, F., Zhang, W., and Hu, B. (2011). A novel approach for determining fatigue resistances of different muscle groups in static cases. Int. J. Ind. Ergon. 41, 10–18. doi: 10.1016/j.ergon.2010.11.005

CrossRef Full Text | Google Scholar

Ma, R., Chablat, D., Bennis, F., Ma, L., Lenarcic, J., and Husty, M. (Eds.). (2012). “Latest advances in robot kinematics,” in Chapter Human Muscle Fatigue Model in Dynamic Motions (Dordrecht: Springer), 349–356. doi: 10.1007/978-94-007-4620-6_44

CrossRef Full Text | Google Scholar

MacIntosh, B. R., Holash, R. J., and Renaud, J. M. (2012). Skeletal muscle fatigue-regulation of excitation-contraction coupling to avoid metabolic catastrophe. J. Cell Sci. 125, 2105–2114. doi: 10.1242/jcs.093674

PubMed Abstract | CrossRef Full Text | Google Scholar

MacIntosh, B. R., and Rassier, D. E. (2002). What is fatigue? Can. J. Appl. Physiol. 27, 42–55. doi: 10.1139/h02-003

PubMed Abstract | CrossRef Full Text | Google Scholar

MacNaughton, M. B., and MacIntosh, B. R. (2007). Impact of length during repetitive contractions on fatigue in rat skeletal muscle. Eur. J. Physiol. 455, 359–366. doi: 10.1007/s00424-007-0273-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Marion, M. S., Wexler, A. S., and Hull, M. L. (2010). Predicting fatigue during electrically stimulated nonisometric contractions. Muscle Nerve 41, 857–867. doi: 10.1002/mus.21603

CrossRef Full Text | Google Scholar

Maton, B., and Gamet, D. (1989). The fatigability of two agonistic muscles in human isometric voluntary submaximal contraction: an EMG study. II. Motor unit firing rate and recruitment. Eur. J. Appl. Physiol. Occupat. Physiol. 58, 369–374. doi: 10.1007/BF00643511

PubMed Abstract | CrossRef Full Text | Google Scholar

Matthijsse, P. C., Hendrich, K. M. M., Rijnsburger, W. H., Woittiez, R. D., and Huijing, P. A. (1987). Ankle angle effects on endurance time, median frequency and mean power of gastrocnemius emg power spectrum: a comparison between individual and group analysis. Ergonomics 30, 1149–1159. doi: 10.1080/00140138708966004

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenzie, D. K., and Gandevia, S. C. (1987). Influence of muscle length on human inspiratory and limb muscle endurance. Respir. Physiol. 67, 171–182. doi: 10.1016/0034-5687(87)90039-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Millar, N. C., and Homsher, E. (1990). The effect of phosphate and calcium on force generation in glycerinated rabbit skeletal muscle fibers. A steady-state and transient kinetic study. J. Biol. Chem. 265, 20234–20240.

PubMed Abstract | Google Scholar

Millard, M., Uchida, T., Seth, A., and Delp, S. L. (2013). Flexing computational muscle: modeling and simulation of musculotendon dynamics. J. Biomech. Eng. 135:12. doi: 10.1115/1.4023390

PubMed Abstract | CrossRef Full Text | Google Scholar

Mörl, F., Siebert, T., and Haeufle, D. F. B. (2016). Contraction dynamics and function of the muscle-tendon complex depend on the muscle fibre-tendon length ratio: a simulation study. Biomech. Model. Mechanobiol. 15, 245–258. doi: 10.1007/s10237-015-0688-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Mörl, F., Siebert, T., Schmitt, S., Blickhan, R., and Günther, M. (2012). Electro-mechanical delay in Hill-type muscle models. J. Mech. Med. Biol. 12, 85–102. doi: 10.1142/S0219519412500856

CrossRef Full Text | Google Scholar

Mosso, A. (1904). Fatigue. London: Swan Sonnenschein.

Google Scholar

Muretta, J. M., Rohde, J. A., Johnsrud, D. O., Cornea, S., and Thomas, D. D. (2015). Direct real-time detection of the structural and biochemical events in the myosin power stroke. Proc. Natl. Acad. Sci. U.S.A. 112, 14272–14277. doi: 10.1073/pnas.1514859112

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakajima, H., Kunioka, Y., Nakano, K., Shimizu, K., Seto, M., and Ando, T. (1997). Scanning force microscopy of the interaction events between a single molecule of heavy meromyosin and actin. Biochem. Biophys. Res. Commun. 234, 178–182. doi: 10.1006/bbrc.1997.6612

PubMed Abstract | CrossRef Full Text | Google Scholar

Narendra, K. S., and Gallman, P. G. (1966). An iterative method for the identification of nonlinear systems using a Hammerstein model. IEEE Trans. Autom. Control 11, 546–550. doi: 10.1109/TAC.1966.1098387

CrossRef Full Text | Google Scholar

Nelson, C. R., and Fitts, R. H. (2014). Effects of low cell pH and elevated inorganic phosphate on the pCa-force relationship in single muscle fibers at near-physiological temperatures. Am. J. Cell Physiol. 306, C670–C678. doi: 10.1152/ajpcell.00347.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Newsholme, E. A. (1971–72). The regulation of phosphofructokinase in muscle. Cardiology 56, 22–34. doi: 10.1159/000169338.

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, A. V., Agre, J. C., Hanson, P., Harrington, M. S., and Nagle, F. J. (1985). Influence of muscle length and force on endurance and pressor responses to isometric exercise. J. Appl. Physiol. 76, 2561–2569. doi: 10.1152/jappl.1994.76.6.2561

PubMed Abstract | CrossRef Full Text | Google Scholar

Nocella, M., Cecchi, G., and Colombini, B. (2017). Phosphate increase during fatigue affects crossbridge kinetics in intact mouse muscle at physiological temperature. J. Physiol. 595, 4317–4328. doi: 10.1113/JP273672

PubMed Abstract | CrossRef Full Text | Google Scholar

Noda, L. (1973). “8 adenylate kinase,” in Group Transfer Part A: Nucleotidyl Transfer Nucleosidyl Transfer Acyl Transfer Phosphoryl Transfer, Volume 8 of The Enzymes, ed P. D. Boyer (New York, NY: Academic Press), 179–305. doi: 10.1016/S1874-6047(08)60068-2

CrossRef Full Text | Google Scholar

Nosek, T. M., Fender, K. Y., and Godt, R. E. (1987). It is diprotonated inorganic phosphate that depresses force in skinned skeletal muscle fibers. Science 236, 191–193. doi: 10.1126/science.3563496

PubMed Abstract | CrossRef Full Text | Google Scholar

Nussbaum, M. A., Clark, L. L., Lanza, M. A., and Rice, K. M. (2001). Fatigue and endurance limits during intermittent overhead work. Am. Ind. Hyg. Assoc. J. 62, 446–456. doi: 10.1202/0002-8894(2001)062<0446:FAELDI>2.0.CO;2

PubMed Abstract | CrossRef Full Text | Google Scholar

Petrofsky, J. S. (1978). Control of the recruitment and firing frequencies of motor units in electrically stimulated muscles in the cat. Med. Biol. Eng. Comput. 16, 302–308. doi: 10.1007/BF02442432

PubMed Abstract | CrossRef Full Text | Google Scholar

Petrofsky, J. S., and Phillips, C. A. (1980). The effect of elbow angle on the isometric strength and endurance of the elbow flexors in men and women. J. Hum. Ergol. 9, 125–131.

PubMed Abstract | Google Scholar

Phillips, S. K., Wiseman, R. W., Woledge, R. C., and Kushmerick, M. J. (1993). The effects of metabolic fuel on force production and resting inorganic phosphate levels in mouse skeletal muscle. J. Physiol. 462, 135–146. doi: 10.1113/jphysiol.1993.sp019547

CrossRef Full Text | Google Scholar

Piazzesi, G., and Lombardi, V. (1995). A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle. Biophys. J. 68, 1966–1979. doi: 10.1016/S0006-3495(95)80374-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Piazzesi, G., Lucii, L., and Lombardi, V. (2002). The size and the speed of the working stroke of muscle myosin and its dependence on the force. J. Physiol. 545, 145–151. doi: 10.1113/jphysiol.2002.028969

PubMed Abstract | CrossRef Full Text | Google Scholar

Ranatunga, K. W. (2010). Force and power generating mechanism(s) in active muscle as revealed from temperature perturbation studies. J. Physiol. 588, 3657–3670. doi: 10.1113/jphysiol.2010.194001

PubMed Abstract | CrossRef Full Text | Google Scholar

Rashedi, E., and Nussbaum, M. A. (2015a). Mathematical models of localized muscle fatigue: sensitivity analysis and assessment of two occupationally-relevant models. PLoS ONE 10:e0143872. doi: 10.1371/journal.pone.0143872

PubMed Abstract | CrossRef Full Text | Google Scholar

Rashedi, E., and Nussbaum, M. A. (2015b). A review of occupationally-relevant models of localised muscle fatigue. Int. J. Hum. Factors Modell. Simul. 5, 61–80. doi: 10.1504/IJHFMS.2015.068119

CrossRef Full Text | Google Scholar

Rassier, D. E. (2000). The effects of length on fatigue and twitch potentiation in human skeletal muscle. Clin. Physiol. 20, 474–482. doi: 10.1046/j.1365-2281.2000.00283.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Reconditi, M. (2006). Recent improvements in small angle X-ray diffraction for the study of muscle physiology. Rep. Prog. Phys. 69, 2709–2759. doi: 10.1088/0034-4885/69/10/R01

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockenfeller, R. (2016). On the application of mathematical methods in hill-type muscle modeling: stability, sensitivity and optimal control (Ph.D. thesis), Universität Koblenz-Landau, Mainz, Germany.

Google Scholar

Rockenfeller, R., and Günther, M. (2016). Extracting low-velocity concentric and eccentric dynamic muscle properties from isometric contraction experiments. Math. Biosci. 278, 77–93. doi: 10.1016/j.mbs.2016.06.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockenfeller, R., and Günther, M. (2017a). Hill equation and Hatze's muscle activation dynamics complement each other: enhanced pharmacological and physiological interpretability of modelled activity-pCa curves. J. Theor. Biol. 431, 11–24. doi: 10.1016/j.jtbi.2017.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockenfeller, R., and Günther, M. (2017b). How to model a muscle's active force-length relation: a comparative study. Comput. Methods Appl. Mech. Eng. 313, 321–336. doi: 10.1016/j.cma.2016.10.003

CrossRef Full Text | Google Scholar

Rockenfeller, R., and Günther, M. (2018). Inter-filament spacing mediates calcium binding to troponin: a simple geometric-mechanistic model explains the shift of force-length maxima with muscle activation. J. Theor. Biol. 454, 240–252. doi: 10.1016/j.jtbi.2018.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Rockenfeller, R., Günther, M., Schmitt, S., and Götz, T. (2015). Comparative sensitivity analysis of muscle activation dynamics. Comput. Math. Methods Med. 2015:585409. doi: 10.1155/2015/585409

PubMed Abstract | CrossRef Full Text | Google Scholar

Rode, C., Siebert, T., and Blickhan, R. (2009). Titin-induced force enhancement and force depression: a ‘sticky-spring' mechanism in muscle contractions? J. Theor. Biol. 259, 350–360. doi: 10.1016/j.jtbi.2009.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Rohmert, W. (1960). Ermittlung von Erholungspausen für statische Arbeit des Menschen (German text). Int. Z. Angew. Physiol. Einschl. Arbeitsphysiol. 18, 123–164. doi: 10.1007/BF00698869

CrossRef Full Text | Google Scholar

Rosenfeld, E. V. (2012). The interrelation between mechanical characteristics of contracting muscle, cross-bridge internal structure, and the mechanism of chemomechanical energy transduction. Eur. Biophys. J. 41, 733–753. doi: 10.1007/s00249-012-0849-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenfeld, E. V., and Günther, M. (2014). An enhanced model of cross-bridge operation with internal elasticity. Eur. Biophys. J. 43, 131–141. doi: 10.1007/s00249-014-0947-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosing, J., and Slater, E. C. (1972). The value of δG° for the hydrolysis of ATP. Biochim. Biophys. Acta 267, 275–290. doi: 10.1016/0005-2728(72)90116-8

CrossRef Full Text | Google Scholar

Sacco, P., McIntyre, D. B., and Jones, D. A. (1985). Effects of length and stimulation frequency on fatigue of the human tibialis anterior muscle. J. Appl. Physiol. 77, 1148–1154. doi: 10.1152/jappl.1994.77.3.1148

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahlin, K., Harris, R. C., and Hultman, E. (1975). Creatine kinase equilibrium and lactate content compared with muscle pH in tissue samples obtained after isometric exercise. Biochem. J. 152, 173–180. doi: 10.1042/bj1520173

PubMed Abstract | CrossRef Full Text | Google Scholar

Saltelli, A., Chan, K., and Scott, E. M. (2000). Sensitivity Analysis. 1st Edn. John Wiley.

Google Scholar

Schappacher-Tilp, G., Leonard, T., Desch, G., and Herzog, W. (2015). A novel three-filament model of force generation in eccentric contraction of skeletal muscles. PLoS ONE 10:e0117634. doi: 10.1371/journal.pone.0117634

PubMed Abstract | CrossRef Full Text | Google Scholar

Scott, S. H., Brown, I. E., and Loeb, G. E. (1996). Mechanics of feline soleus: I. Effect of fascicle length and velocity on force output. J. Muscle Res. Cell Motil. 17, 207–219. doi: 10.1007/BF00124243

PubMed Abstract | CrossRef Full Text | Google Scholar

Seow, C. Y. (2013). Hill's equation of muscle performance and its hidden insight on molecular mechanisms. J. Gen. Physiol. 142, 561–573. doi: 10.1085/jgp.201311107

PubMed Abstract | CrossRef Full Text | Google Scholar

Seth, D., Chablat, D. C., Bennis, F., Sakka, S., Jubeau, M., and Nordez, A. (2016). “New dynamic muscle fatigue model to limit musculo-skeletal disorder,” in The 2016 Virtual Reality International Conference (New York, NY: Association for Computing Machinery), 9. doi: 10.1145/2927929.2927935

CrossRef Full Text | Google Scholar

Shorten, P. R., O'Callaghan, P., Davidson, J. B., and Soboleva, T. K. (2007). A mathematical model of fatigue in skeletal muscle force contraction. J. Muscle Res. Cell Motil. 28, 293–313. doi: 10.1007/s10974-007-9125-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Siebert, T., Leichsenring, K., Rode, C., Wick, C., Stutzig, N., Schubert, H., et al. (2015). Three-dimensional muscle architecture and comprehensive dynamic properties of rabbit gastrocnemius, plantaris and soleus: input for simulation studies. PLoS ONE 10:e0130985. doi: 10.1371/journal.pone.0130985

PubMed Abstract | CrossRef Full Text | Google Scholar

Siebert, T., Rode, C., Herzog, W., Till, O., and Blickhan, R. (2008). Nonlinearities make a difference: comparison of two common Hill-type models with real muscle. Biol. Cybern. 98, 133–143. doi: 10.1007/s00422-007-0197-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Siebert, T., Sust, M., Thaller, S., Tilp, M., and Wagner, H. (2007). An improved method to determine neuromuscular properties using force laws–from single muscle to applications in human movements. Hum. Mov. Sci. 26, 320–341. doi: 10.1016/j.humov.2007.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Sieck, G. C., Ferreira, L. F., Reid, M. V., and Mantilla, C. B. (2013). Mechanical properties of respiratory muscles. Compr. Physiol. 3, 1553–1567. doi: 10.1002/cphy.c130003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sih, B., Ng, L., and Stuhmiller, J. (2012). Generalization of a model based on biophysical concepts of muscle activation, fatigue and recovery that explains exercise performance. Int. J. Sports Med. 33, 258–267. doi: 10.1055/s-0031-1297958

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, L. A., Chock, P. B., and Eisenberg, E. (1981). Mechanism of the actomyosin ATPase: effect of actin on the ATP hydrolysis step. Proc. Natl. Acad. Sci. U.S.A. 78, 1346–1350. doi: 10.1073/pnas.78.3.1346

PubMed Abstract | CrossRef Full Text | Google Scholar

Stutzig, N., Rzanny, R., Moll, K., Gussew, A., Reichenbach, J. R., and Siebert, T. (2017). The pH heterogeneity in human calf muscle during neuromuscular electrical stimulation. Magn. Reson. Med. 77, 2097–2106. doi: 10.1002/mrm.26329

PubMed Abstract | CrossRef Full Text | Google Scholar

Suzuki, Y., Yasunaga, T., Ohkura, R., Wakabayashi, T., and Sutoh, K. (1998). Swing of the lever arm of a myosin motor at the isomerization and phosphate-release steps. Nature 396, 380–383. doi: 10.1038/24640

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, Y., Shuman, H., and Goldman, Y. E. (2004). Coupling between phosphate release and force generation in muscle actomyosin. Philos. Trans. R. Soc. B Biol. Sci. 359, 1913–1920. doi: 10.1098/rstb.2004.1561

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, C. Y., Stojanovic, B., Tsui, C. P., and Kojic, M. (2005). Modeling of muscle fatigue using Hill's model. Biomed. Mater. Eng. 15, 341–348.

PubMed Abstract | Google Scholar

Tesch, P., Sjödin, B., Thorstensson, A., and Karlsson, J. (1978). Muscle fatigue and its relation to lactate accumulation and LDH activity in man. Acta Physiol. Scand. 103, 413–420. doi: 10.1111/j.1748-1716.1978.tb06235.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tesi, C., Barman, T., and Lionne, C. (2017). Are there two different binding sites for atp on the myosin head, or only one that switches between two conformers? J. Muscle Res. Cell Motil. 38, 137–142. doi: 10.1007/s10974-017-9480-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tesi, C., Colomo, F., Piroddi, N., and Poggesi, C. (2002). Characterization of the cross-bridge force-generating step using inorganic phosphate and BDM in myofibrils from rabbit skeletal muscles. J. Physiol. 541, 187–199. doi: 10.1113/jphysiol.2001.013418

PubMed Abstract | CrossRef Full Text | Google Scholar

Till, O., Siebert, T., Rode, C., and Blickhan, R. (2008). Characterization of isovelocity extension of activated muscle: a hill-type model for eccentric contractions and a method for parameter determination. J. Theor. Biol. 225, 176–187. doi: 10.1016/j.jtbi.2008.08.009

CrossRef Full Text | Google Scholar

Tonomura, Y., Nakamura, H., Kinoshita, N., Onishi, H., and Shigekawa, M. (1966). The pre-steady state of the myosin-adenosine triphosphate system. J. Biochem. 66, 599–618.

PubMed Abstract | Google Scholar

Trentham, D. R., Eccleston, J. F., and Bagshaw, C. R. (1976). Kinetic analysis of ATPase mechanisms. Q. Rev. Biophys. 9, 217–281. doi: 10.1017/S0033583500002419

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dieën, J. H., and oude Vrielink, H. H. E. (1994). The use of the relation between relative force and endurance time. Ergonomics 37, 231–243. doi: 10.1080/00140139408963641

PubMed Abstract | CrossRef Full Text | Google Scholar

Vøllestad, N. K. (1997). Measurement of human muscle fatigue. J. Neurosci. Methods 74, 219–227. doi: 10.1016/S0165-0270(97)02251-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner, H., Siebert, T., Ellerby, D. J., Marsh, R. L., and Blickhan, R. (2005). ISOFIT: a model-based method to measure muscle-tendon properties simultaneously. Biomech. Model. Mechanobiol. 4, 10–19. doi: 10.1007/s10237-005-0068-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L. C., and Kernell, D. (2001). Fibre type regionalisation in lower hindlimb muscles of rabbit, rat and mouse: a comparative study. J. Anat. 199, 631–643. doi: 10.1046/j.1469-7580.2001.19960631.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Westerblad, H., and Allen, D. G. (1991). Changes of myoplasmic calcium concentration during fatigue in single mouse muscle fibers. J. Gen. Physiol. 98, 615–635. doi: 10.1085/jgp.98.3.615

PubMed Abstract | CrossRef Full Text | Google Scholar

Westerblad, H., and Allen, D. G. (1993). The contribution of [Ca²⁺]i to the slowing of relaxation in fatigued single fibres from mouse skeletal muscle. J. Physiol. 468, 729–740. doi: 10.1113/jphysiol.1993.sp019797

CrossRef Full Text | Google Scholar

Westerblad, H., and Allen, D. G. (2002). Recent advances in the understanding of skeletal muscle fatigue. Curr. Opin. Rheumatol. 14, 648–652. doi: 10.1097/00002281-200211000-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Westerblad, H., Allen, D. G., and Lännergren, J. (2002). Muscle fatigue: lactic acid or inorganic phosphate the major cause? News Physiol. Sci. 17, 17–21. doi: 10.1152/physiologyonline.2002.17.1.17

PubMed Abstract | CrossRef Full Text | Google Scholar

Westerblad, H., Bruton, J. D., and Lännergren, J. (1997). The effect of intracellular pH on contractile function of intact, single fibres of mouse muscle declines with increasing temperature. J. Physiol. 500, 193–204. doi: 10.1113/jphysiol.1997.sp022009

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiener, N. (1942). Response of a Nonlinear System to Noise. Restricted report v-16, no 129, 112pp. Declassified July 1946, published as rep. no. pb-1-58087, Technical report, Radiation Lab MIT.

Google Scholar

Willems, M. E. T., and Stauber, W. T. (2002). Fatigue and recovery at long and short muscle lengths after eccentric training. Med. Sci. Sports Exerc. 34, 1738–1743. doi: 10.1097/00005768-200211000-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wills, A., Schön, T. B., Ljung, L., and Ninness, B. (2013). Identification of hammerstein-wiener models. Automatica 49, 70–81. doi: 10.1016/j.automatica.2012.09.018

CrossRef Full Text | Google Scholar

Xia, T., and Frey-Law, L. A. (2008). A theoretical approach for modeling peripheral muscle fatigue and recovery. J. Biomech. 41, 3046–3052. doi: 10.1016/j.jbiomech.2008.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X. C., and Feng, W. (2016). Thermodynamic aspects of ATP hydrolysis of actomyosin complex. Biophys. Rev. 2, 87–94. doi: 10.1007/s41048-016-0032-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendices

A. Initial muscle model

A.0.1. Model nomenclature

The herein sketched macroscopic Hill-type muscle model is completely based on prior publications (Günther et al., 2007; Haeufle et al., 2014b; Mörl et al., 2012), see also (Rockenfeller and Günther, 2016, App. A). Let ℓ_X and F_X denote the length and force, respectively, of element X, which represents either the whole muscle tendon unit (MTU), the fiber material (contractile element, CE), the connective tissue (parallel elastic element, PEE), or the tendon (serial elastic element, SEE, and serial damping element, SDE).

A.0.2. Model assumptions

It holds that ℓ_CE = ℓ_PEE, as well as ℓ_SEE = ℓ_SDE, and ℓ_CE+ℓ_SEE = ℓ_MTU. Further, the force equilibrium

$\begin{array}{l}{F}_{\text{MTU}}=F_{\text{CE}}\left(\tilde{q},{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}\right)+F_{\text{PEE}}({\ell }_{\text{CE}})=F_{\text{SEE}}({\ell }_{\text{CE}},{\ell }_{\text{MTU}})\\ +F_{\text{SDE}}\left(\tilde{q},{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{MTU}}\right)& (9)\end{array}$

is assumed to hold at any time instance, see (Haeufle et al., 2014b, Figure 1) for a model sketch.

A.0.3. Activation dynamics

After a piece-wise continuous neural impulse u = u(t) at the neuromuscular junction, the relative calcium ion concentration $\tilde{c}$ within the sarcolemma changes with a time constant m via the ODE

$\begin{array}{ll}\frac{\text{d}}{\text{d}t}\tilde{c}=m·(u-\tilde{c}).& (10)\end{array}$

Subsequently, calcium binds to troponin C and causes tropomyosin to move aside from the actin filament, thereby uncovering the seven active sites (Reconditi, 2006). The relative amount of calcium-bound troponin is commonly referred to as (troponin-)activity (Hatze, 1977, 1978, 1981) and denoted $\tilde{q}$. This process, herein referred to by the term activation dynamics, is known to be length-dependent (Kistemaker et al., 2005):

$\begin{array}{ll}\tilde{q}(\tilde{c},{\tilde{\ell }}_{\text{CE}})=\frac{{\tilde{q}}_{min}+{\left(\varpi ({\tilde{\ell }}_{\text{CE}})·\tilde{c}\right)}^{\nu }}{1+{\left(\varpi ({\tilde{\ell }}_{\text{CE}})·\tilde{c}\right)}^{\nu }},\text{ where }\varpi ({\tilde{\ell }}_{\text{CE}})={\varpi }_{\text{opt}}·{\tilde{\ell }}_{\text{CE}}.& (11)\end{array}$

The formulation includes a minimum troponin-activity ${\tilde{q}}_{min}$ of the muscle at rest and a linear length dependency function ϖ describing the troponin C density with varying interfilamentary spacing (Rockenfeller and Günther, 2018), based on sarcomere volume constancy (Dragomir, 1970). The Hatze exponent ν has previously (Rockenfeller and Günther, 2017a) been equated to the exponent of the Hill equation (Hill, 1910). However, so far no mechanistic explanation has been provided on how to convert the relative amount of calcium-bound troponin (troponin-activity: $\tilde{q}$) to the relative amount of built-up cross-bridges (actual activity: ã), i.e., it is implied that myosin heads will automatically bind to available active sites and proceed with force generation. In section 3, we investigated this missing link of cross-bridge formation from a physiological point of view and incorporate the concept of phosphate dynamics. For this summary, the symbol ã is used for overall activity and the reader may keep in mind that it incorporates both, activation and phosphate dynamics.

A.0.4. Force-length-velocity relation

Depending on filamentary overlap, the muscle fibers are assumed to produce a force proportional to the isometric force-length relation (FLR), cf. (Günther et al., 2007; Rockenfeller and Günther, 2017b):

$\begin{array}{ll}{\tilde{F}}_{\text{isom}}({\tilde{\ell }}_{\text{CE}})=\{\begin{array}{ll}\exp \left(-{\left|\frac{{\tilde{\ell }}_{\text{CE}}-1}{\Delta W_{\text{asc}}}\right|}^{v_{\text{asc}}}\right),\hfill & \text{if }{\tilde{\ell }}_{\text{CE}}\le 1\hfill \\ \exp \left(-{\left|\frac{{\tilde{\ell }}_{\text{CE}}-1}{\Delta W_{\text{des}}}\right|}^{v_{\text{des}}}\right),\hfill & \text{if }{\tilde{\ell }}_{\text{CE}}>1\hfill \end{array}.& (12)\end{array}$

The parameters ΔW_asc, ΔW_des, and ν_asc, ν_des represent the width and kurtosis of the double-exponential function, respectively. The optimal fiber length ℓ_{CE, opt}, at which maximum force F_max can be produced, represents the situation in the fully activated muscle tissue. It is known that in sub-maximally activated fibers the force maximum shifts to longer length with decreasing activity (Rockenfeller and Günther, 2018). The resulting, velocity-dependent fiber force is modeled according to the Hill relation (Hill, 1938):

$\begin{array}{l}{F}_{\text{CE}}\left(\tilde{a},{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}\right)\\ =\{\begin{array}{ll}{F}_{\max }·\left(\frac{\tilde{a}·{\tilde{F}}_{\text{isom}}+A_{\text{rel}}}{1-\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}/(B_{\text{rel}}\cdot {\ell }_{\text{CE},\text{opt}})}-A_{\text{rel}}\right),\hfill & \text{if}\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}\le 0\hfill \\ F_{\max }·\left(\frac{\tilde{a}·{\tilde{F}}_{\text{isom}}+A_{\text{rel,e}}}{1-\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}/(B_{\text{rel,e}}\cdot {\ell }_{\text{CE},\text{opt}})}-A_{\text{rel,e}}\right),\hfill & \text{if}\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}>0\hfill \end{array},& (13)\end{array}$

where in the concentric case the Hill parameters A_rel, B_rel are activity- and length-dependent

$\begin{array}{ll}{A}_{\text{rel}}(\tilde{a},{\ell }_{\text{CE}})=\{\begin{array}{ll}{A}_{\text{rel},0}\cdot (1+3\cdot \tilde{a})/4,\hfill & \text{ if }{\tilde{\ell }}_{\text{CE}}\le 1\hfill \\ A_{\text{rel},0}\cdot {\tilde{F}}_{\text{isom}}\cdot (1+3\cdot \tilde{a})/4,\hfill & \text{ if }{\tilde{\ell }}_{\text{CE}}>1\hfill \end{array},& (14)\end{array}$

$\begin{array}{ll}{B}_{\text{rel}}(\tilde{a})=B_{\text{rel,0}}·(3+4·\tilde{a})/7.& (15)\end{array}$

In the eccentric case, the Hill hyperbola is reflected to account for experimental findings (Katz, 1939; Komi, 1973; Till et al., 2008):

$\begin{array}{ll}{A}_{\text{rel,e}}(\tilde{a},{\ell }_{\text{CE}})=-F_{\text{e}}·\tilde{a}·{\tilde{F}}_{\text{isom}},B_{\text{rel,e}}(\tilde{a},{\ell }_{\text{CE}})=\frac{B_{\text{rel}}·(1-F_{\text{e}})}{S_{\text{e}}·\left(1+\frac{A_{\text{rel}}}{\tilde{a}·{\tilde{F}}_{\text{isom}}}\right)}.& (16)\end{array}$

The shape parameters F_e, S_e represent the relative eccentric force limit and the degree of non-differentiability (in terms of a Lipschitz constant) at $\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0$ (Haeufle et al., 2014a, Equation 9), respectively.

A.0.5. Passive components:

Parallel and serial elastic force-length relations are considered as non-linear and nonlinear-linear elastic springs, respectively:

$\begin{array}{l}{F}_{\text{PEE}}({\tilde{\ell }}_{\text{CE}})\\ =\{\begin{array}{ll}0,\hfill & \text{ if }{\tilde{\ell }}_{\text{CE}}\le L_{\text{PEE},0}\hfill \\ F_{\text{PEE}}\cdot F_{\max }\cdot {\left(\frac{{\tilde{\ell }}_{\text{CE}}-L_{\text{PEE,0}}}{\text{Δ}{W}_{\text{des}}+1-L_{\text{PEE,0}}}\right)}^{{}^v\text{P}\text{EE}},\hfill & \text{ if }{\tilde{\ell }}_{\text{CE}}>L_{\text{PEE},0}\hfill \end{array},& (17)\end{array}$

$\begin{array}{l}{F}_{\text{SEE}}({\ell }_{\text{CE}},{\ell }_{\text{MTU}})\\ =\{\begin{array}{ll}0,\hfill & \text{ if }{\tilde{\ell }}_{\text{SEE}}\le 1\hfill \\ \text{Δ}{F}_{\text{SEE},0}\cdot {\left(\frac{{\tilde{\ell }}_{\text{SEE}}-1}{\text{Δ}{U}_{\text{SEE},\text{nll}}}\right)}^{\text{Δ}{U}_{\text{SEE},\text{nll}}/\text{Δ}{U}_{\text{SEE},\text{l}}},\hfill & \text{ if }{\tilde{\ell }}_{\text{SEE}}\le 1+\text{Δ}{U}_{\text{SEE},\text{nll}}\hfill \\ \text{Δ}{F}_{\text{SEE},0}\cdot \left(1+\frac{{\tilde{\ell }}_{\text{SEE}}-1+\text{Δ}{U}_{\text{SEE},\text{nll}}}{\text{Δ}{U}_{\text{SEE},\text{l}}}\right),\hfill & \text{ if }{\tilde{\ell }}_{\text{SEE}}>1+\text{Δ}{U}_{\text{SEE},\text{nll}}\hfill \end{array}.& (18)\end{array}$

The PEE slack length as well as characteristic force is given relative to the optimal fiber length and maximum fiber force (${{L}}_{\text{PEE,0}}$ and ${{F}}_{\text{PEE}}$), whereas the SEE slack length ℓ_{SEE, 0} and force shape parameter ΔF_{SEE, 0} are considered as parameters. The relative SEE length ${\tilde{\ell }}_{\text{SEE}}$ is obtained by ${\tilde{\ell }}_{\text{SEE}}:=({\ell }_{\text{MTU}}-{\ell }_{\text{CE}})/{\ell }_{\text{SEE,0}}$. The slopes of the nonlinear and linear part of the serial elastic FLR are represented by ΔU_{SEE, nll} and ΔU_{SEE, l}, respectively, see also (Günther et al., 2007, Figure 4). The remaining serial damping element is considered to be linearly dependent on MTU force,

$\begin{array}{l}{F}_{\text{SDE}}\left(\tilde{a},{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}},\frac{\text{d}}{\text{d}t}{\ell }_{\text{MTU}}\right)=D_{\text{SE}}·\frac{F_{max}·A_{\text{rel,0}}}{B_{\text{rel,0}}·{\ell }_{\text{CE,opt}}}·\\ \left((1-R_{\text{SE}})·\frac{F_{\text{CE}}+F_{\text{PEE}}}{F_{max}}+R_{\text{SE}}\right)·\left(\frac{\text{d}}{\text{d}t}{\ell }_{\text{MTU}}-\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}\right),& (19)\end{array}$

where D_SE represents the slope and R_SE the minimum damping coefficient.

A.0.6. Contraction dynamics

Solving the force equilibrium (9) with respect to $\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}$ yields a quadratic equation, which can be transferred to a nonlinear ODE

$\begin{array}{ll}\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=f(\tilde{a},{\ell }_{\text{CE}})& (20)\end{array}$

that is referred to as the term contraction dynamics. For details see (Haeufle et al., 2014b, sections 2.5 and 2.6). Note that in the stationary case ($\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0$) the MTU force from Equation (9) is given by

$\begin{array}{ll}{F}_{\text{MTU}}(\tilde{a},{\ell }_{\text{CE}}){\mid }_{\frac{\text{d}}{\text{d}t}{\ell }_{\text{CE}}=0}=\tilde{a}·F_{max}·{\tilde{F}}_{\text{isom}}({\ell }_{\text{CE}})+F_{\text{PEE}}({\ell }_{\text{CE}}).& (21)\end{array}$

B. ATP reaction kinetics

Basis for our kinetic model of phosphate dynamics is the following reaction scheme (Allen and Orchard, 1987)

$\begin{array}{ll}\text{ATP}<=>\text{ADP}+{\text{P}}_{\text{i}},& (22)\end{array}$

$\begin{array}{ll}\text{PCr}+\text{ADP}<=>\text{ATP + Cr},& (23)\end{array}$

$\begin{array}{ll}2\text{ADP}<=>\text{ATP + AMP},& (24)\end{array}$

including the auxiliary constraints

$\begin{array}{l}\text{total adenine concentration:}{c}_{\text{ad}}=7\text{mM}=[\text{ATP}]+[\text{ADP}]\\ +[\text{AMP}],& (25)\end{array}$

$\begin{array}{ll}\text{total creatine concentration:}{c}_{\text{cr}}=25\text{mM}=[\text{PCr}]+[Cr],& (26)\end{array}$

$\begin{array}{l}\text{total phosphate concentration:}{c}_{\text{ph}}=46\text{mM}=3[\text{ATP}]\\ +2[\text{ADP}]+[\text{AMP}]\\ +[\text{PCr}]+[{\text{P}}_i].& (27)\end{array}$

For the equilibrium reactions (23) and (24) the law of mass action states

$\begin{array}{ll}{K}_{\text{ATP}}=\frac{\text{[ATP]}\cdot c_0}{\text{[ADP]}\cdot [{\text{P}}_i]}=\frac{k_{\text{con}}}{k_{\text{hyd}}},& (28)\end{array}$

$\begin{array}{ll}{K}_{\text{crk}}=\frac{[\text{ATP}][\text{Cr}]}{[\text{ADP}][\text{PCr}]}=200,& (29)\end{array}$

$\begin{array}{ll}{K}_{\text{adk}}=\frac{[\text{ATP}][\text{AMP}]}{{[\text{ADP}]}^2}=1.& (30)\end{array}$

where k_hyd and k_con denote the ATP hydrolysis and condensation rate constant, respectively. The equilibrium constants for creatine kinase (crk) and adenylate kinase (adk) are likewise taken from Allen and Orchard (1987). It should be noted that there exist other sources with different values e.g., for initial [PCr] (Allen et al., 2008; Karatzaferi et al., 2001), as well as the equilibrium constants for the creatine kinase (Bechtel and Best, 1985; Gadian et al., 1981; Lawson and Veech, 1979; Sahlin et al., 1975) and the adenylate kinase (Kübler and Katz, 1977; Linari et al., 2010; Newsholme, 1 72; Noda, 1973) reaction. For the sake of consistency, however, we obtained all values from the same source, because the presented data do not allow a detailed assessment and thus estimation of these parameters.

For the pseudo-first order (i.e., the concentration of water is assumed to remain constant) reaction [Equations (22) or (1) together with (28)] it holds that

$\begin{array}{ll}(-\frac{\text{d}}{\text{d}t}[\text{ATP}](t)=\frac{\text{d}}{\text{d}t}[\text{ADP}](t)=\frac{\text{d}}{\text{d}t}[P_i](t)& (31)\end{array}$

$\begin{array}{ll}=k_{\text{hyd}}\left(q,{\tilde{\ell }}_{\text{CE}}\right)\cdot [\text{ATP}]\cdot c_0-k_{\text{con}}\cdot [\text{ADP}]\cdot [{\text{P}}_{i\prime }]& (32)\end{array}$

where [ATP](t) and [P_i](t) denote the time evolution of [ATP] and [P_i]. Note that k_hyd is assumed to depend on the current amount of calcium-bound troponin q and the current CE length ${\tilde{\ell }}_{\text{CE}}$, cf. section 3. As an initial condition, set [ATP](0) = [ATP]₀. As the equilibrium constitutes for a rather lumped description of the Lymn-Taylor cycle (see Discussion 5.2), the choice of a suitable value for k_con constitutes for the limiting factor of the reaction, approximately 0.1 Hz (Linari et al., 2010).

In order to obtain a numerical solution for [ATP](t), plug Equation (25) in Equations (27) and (30) as well as Equation (26) in Equation (29), yielding

$\begin{array}{ll}[\text{ADP}]=\frac{-[\text{ATP}]+\sqrt{4·K_{\text{adk}}·c_{\text{ad}}·[\text{ATP}]-(4·K_{\text{adk}}-1)·{[\text{ATP}]}^2}}{2·K_{\text{adk}}},& (33)\end{array}$

$\begin{array}{ll}[\text{AMP}]=c_{\text{ad}}-[\text{ATP}]-[\text{ADP}],& (34)\end{array}$

$\begin{array}{ll}[\text{PCr}]=\frac{c_{\text{cr}}·[\text{ATP}]}{[\text{ATP}]+K_{\text{crk}}·[\text{ADP}]},& (35)\end{array}$

$\begin{array}{ll}[\text{Cr}]=c_{\text{cr}}-[\text{PCr}].& (36)\end{array}$

The hereof obtained expressions for [PCr] (Equation 35) and [ADP] (Equation 33) are plugged in Equation (27), resulting in variable [ATP]:

$\begin{array}{ll}0=c_{\text{ad}}-c_{\text{ph}}+[{\text{P}}_i]+2\cdot [\text{ATP}]+[\text{ADP}]+[\text{PCr}],& (37)\end{array}$

thus expressing [ATP] as a function of the only remaining variable [P_i]. In fact, both Equations (32) and (37) imply the same property of the reaction scheme (22)-(30), namely, that [ATP], [ATP] or [P_i] could either be treated as a new state variable: if one is known, the others can be calculated. Here, the root [ATP]([P_i]) of Equation (37), with additionally using Equations (33) and (35), feeds the right side of Equation (32), which enables to continuously update [P_i] by integrating the ODE. For these calculations, self-implemented Runge-Kutta and bisection methods were used. Initial conditions of all reactants followed from assuming [ADP]₀ = 0.01 mM (Allen et al., 2008) in the resting fiber: [ATP]₀ = 6.99 mM, ${\left[AMP\right]}_0=10^{-5}\text{mM}$, [PCr]₀ = 19.45 mM, [Cr]₀ = 5.55 mM , [P_i]₀ = [Cr]₀ + [ADP]₀ = 5.56 mM. Further from physiological experiments, a minimum of remaining ATP, namely [ATP]_min = 1.2 mM was found in exhausted fibers (Allen and Orchard, 1987) and included as a constraint in the bisection method. Resulting time courses of most involved reactants, including their length-dependency, are displayed in Figure A1 in Supplementary Material.

The addressed interaction between ATP, ADP, and P_i is known to influence the chemical potential (Zhang and Feng, 2016) or affinity (Allen and Orchard, 1987; Cooke, 2007; Hancock et al., 2005) of ATP, write μ_ATP, which constitutes for the new state variable of phosphate dynamics. In a nutshell, the position of equilibrium (22) influences the amount of free Gibbs energy from ATP hydrolysis via

$\begin{array}{ll}{\mu }_{\text{ATP}}:=-\text{Δ}{G}_{\text{ATP}}^{°}+R\cdot T\cdot \ln \left(\frac{[\text{ATP}]·c_0}{[\text{ADP}]\cdot [{\text{P}}_i]}\right),& (38)\end{array}$

where $\Delta G_{\text{ATP}}^{°}\approx -30$ kJ/mol (Allen and Orchard, 1987; Cooke, 2007; Guynn and Veech, 1973; Rosing and Slater, 1972; Zhang and Feng, 2016) denotes the standard free enthalpy, R = 8.31·10⁻³ kJ/(mol·K) the ideal gas constant, T = 303 K a temperature of 30°C, and c₀ = 1 M the standard concentration. The latter serves as a normalization factor. An estimate for the theoretically maximum value for the chemical potential in the resting fiber can be given by

$\begin{array}{ll}{\mu }_{\text{ATP},\text{max}}=-\text{Δ}{G}_{\text{ATP}}^{°}+R\cdot T\cdot \ln \left(\frac{{[\text{ATP}]}_0\cdot c_0}{{[\text{ADP}]}_0\cdot {[P_i]}_0}\right)\approx 59.6\frac{\text{kJ}}{\text{mol}}& (39)\end{array}$

in accordance with literature data (Barclay, 2015). The relative (normalized) ATP affinity is then defined as ${\tilde{\mu }}_{\text{ATP}}:={\mu }_{\text{ATP}}/{\mu }_{\text{ATP,max}}$, for which a time course is likewise displayed in Figure A1 in Supplementary Material.

The time derivative of ${\tilde{\mu }}_{\text{ATP}}$ as used in Equation (8) can finally be calculated as

$\begin{array}{ll}\frac{\text{d}}{\text{d}t}{\tilde{\mu }}_{\text{ATP}}=\frac{R\cdot T}{{\mu }_{\text{ATP},\text{max}}}\cdot \left(\frac{\frac{\text{d}}{\text{d}t}[\text{ATP}]}{[\text{ATP}]}-\frac{\frac{\text{d}}{\text{d}t}[\text{ADP}]}{[\text{ADP}]}-\frac{\frac{\text{d}}{\text{d}t}[P_i]}{[{\text{P}}_i]}\right)~-k_{\text{hyd}},& (40)\end{array}$

where the latter proportionality holds after inserting Equation (32) in Equation (40) and assuming a constant concentration of the involved phosphate compounds.

C. Parameter estimation and sensitivity analysis

For the mathematical formalization of the introduced model, define the domain of the neural stimulation u = u(t) as ${U}:=C^{pw}(\left[t_{start},t_{end}\right],\left[0,1\right])$, i.e., the set of piece-wise continuous functions mapping the experimentally covered time horizon [t_start, t_end] (here: t_start = −0.1s and t_end = 1.1s) to the interval [0, 1]. Let ${\Lambda }_M\in {\text{R}}^d$ denote the tuple of model parameters, representing the corresponding muscle M∈{GAS, PLA}. Denote the modeled time evolution of muscle force by

$\begin{array}{ll}{F}_M:{U}\times {\text{R}}^d\to C([t_{start},t_{end}],\text{R}),F_M(u,{\Lambda }_M)=F_M(t).& (41)\end{array}$

Our aim is to find the optimal parameter set ${\Lambda }_M^{*}$ such that for a given experimental stimulation protocol (input) $\overline{u}$ the measured force (data) ${\overline{F}}_M(t)$ is at best approximated by the model force (output) F_M(t) in a least-squares sense:

$\begin{array}{ll}{\Lambda }_M^{*}={\displaystyle \underset{{\Lambda }_M\in {\text{R}}^d}{arg min}}{J}({\Lambda }_M):=||F_M(\overline{u},{\Lambda }_M)-{\overline{F}}_M(t){||}_2.& (42)\end{array}$

For solving problem (42), we used a pre-implemented MatLab (R2018b, The MathWorks, Natick, USA) routine, namely lsqcurvefit, utilizing a trust-region-reflective algorithm.

In order to asses the influence of every model parameter at each time instance, additional sensitivity values are provided in section 4.2. The local sensitivity (Rockenfeller et al., 2015; Saltelli et al., 2000) of a certain parameter λ ∈ Λ_M is thereby defined as the partial derivative of the objective function with respect to λ at the optimum, i.e., $S_{M,\lambda }(t):=\partial {J}({\Lambda }_M^{*})/\partial \lambda$. A large positive (negative) sensitivity means that only relatively small changes in λ to larger (smaller) values are needed in order to attune model output and measured data. The smaller the sensitivity the larger the required relative change. Hence, parameters with high absolute sensitivities have to be treated with high accuracy, but are easier determinable than insensitive parameters. Conversely, the model is not prone to uncertainties in parameters with comparatively low absolute sensitivities, but these parameters are hardly determinable, cf. (Rockenfeller, 2016, Ch. 5.3) and (Rockenfeller and Günther, 2016).

FIGURE A1

Figure A1. Time courses of [ATP], [P_i], [ADP], and [PCr] (left ordinate) as well as ${\tilde{\mu }}_{\text{ATP}}$ (right ordinate) for GAS (A) and PLA (B). Solid lines represent the time courses for the shortest muscle examined in our experiments (cf. Figure 1). Dashed lines represent the respective longest muscle. Time courses for [Cr] are not shown due to a near congruence with [P_i](t). Note that the ATP affinity decreases substantially, due to increasing P_i, while the [ATP] itself remains almost constant.

List of symbols and abbreviations