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) 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 dynamics (or short phosphate dynamics) (Allen and Orchard, 1987) and the relative chemical potential of ATP as the corresponding state variable that is diminished by both ATP depletion as well as P 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, Fmax) to its plateau region (100% Fmax). 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

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).

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.

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).

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

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 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] 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

namely ATP hydrolysis/condensation. Herein, kcon and khyd 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

where c0 = 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] can be expressed in terms of the first order ODE

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.

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 and M), 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.

For a start, in our model, KATP and in particular khyd, are assumed to be dependent on the current amount of calcium-bound troponin C terminals as well as relative fiber length , cf. Appendix A. The condensation rate kcon is assumed to be a constant. The dependence of khyd on 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 representing minimum troponin-activity. Hence, the rate constant khyd is assumed to increase along with , 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 release (Bickham et al., 2011, Figure 3a) are inversely proportional to fiber length (see also section 5.5). Summarizing, the following ansatz is obtained:

with the parameter 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

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

a rising level of P (and simultaneously of ADP) deteriorates the relative chemical potential of ATP, denoted , which is defined as

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 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.,

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 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

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

Symbol	Value	Unit	Source	Meaning
[ADP]0	0.01	[mM]	Allen et al., 2008	Initial (t = 0) [ADP]
[ATP]min	1.2	[mM]	Allen et al., 2008	Minimum [ATP]
cad	7	[mM]	Allen and Orchard, 1987	Total adenine concentration
ccr	25	[mM]	Allen and Orchard, 1987	Total creatine concentration
cph	46	[mM]	Allen and Orchard, 1987	Total phosphate concentration
Kadk	1	[ ]	Allen and Orchard, 1987	Adenylate kinase equilibrium constant
Kck	200	[]	Allen and Orchard, 1987	Creatine kinase equilibrium constant
kcon	0.1	[Hz]	Linari et al., 2010	ATP condensation rate constant
k^_hyd	Open	[Hz]	–	ATP hydrolysis rate constant

All values for the eight fixed parameters were taken from literature. The optimized values for the open parameter .

4. Results

In total, nGAS = 7 and nPLA = 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

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., ).

Table 2

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

Parameter			Model		Parameter	Model
Symbol	Unit	ΛGASbounds	GAS−	GAS+	ΛPLAbounds	PLA−	PLA+
ℓSEE, 0	[m]	[0.09…0.11]	0.0996	0.1017	[0.09…0.11]	0.1045	0.1054
ΔUSEE, nll	[]	[0.04…0.11]	0.103	0.0785	[0.04…0.11]	0.0620	0.0762
ΔFSEE, 0	[]	[80…160]	123	109	[50…70]	69.4	68.1
ΔUSEE, l	[]	[0.01…0.08]	0.0524	0.0482	[0.01…0.08]	0.0409	0.0414
F_PEE	[]	[0…1.2]	0.133	0.135	[0…1.2]	0.118	0.0908
L_PEE	[]	[0.9…1.5]	1.19	1.02	[1…1.5]	1.00	1.00
νPEE	[]	[2.5…5]	5.00	3.16	[2.5…5]	2.50	2.50
ΔWasc	[]	[0.1…1]	0.215	0.404	[0.1…1]	0.237	0.268
ΔWdes	[]	[0.1…1]	0.690	0.330	[0.1 …1]	0.429	0.570
νasc	[]	[1.5…8]	4.10	6.75	[1.5…5]	1.57	1.84
νdes	[]	[1.5…6]	4.46	2.36	[1.5…5]	3.89	2.83
Fmax	[N]	[80…160]	120	133	[50…70]	57.1	63.7
ℓCE, opt	[m]	[0.014…0.025]	0.0160	0.0198	[0.005…0.02]	0.0130	0.0106
Arel, 0	[]	[0.03…0.5]	0.0786	0.0655	[0.03…0.5]	0.244	0.240
Brel, 0	[1/s]	[1…10]	6.30	2.91	[1…10]	10.0	6.58
DSDE	[]	[0.1…10]	1.65	0.149	[0.1…10]	10.0	3.88
RSDE	[]	[0.01…1]	0.0233	0.0100	[0.1…1]	1.00	0.356
Se	[]	[1…4.5]	1.00	2.85	[1…4.5]	4.50	3.40
Fe	[]	[1.1…2]	1.10	1.14	[1.1…2]	1.61	1.32
q~_min	[]	[0.001…0.01]	0.00576	0.00944	[0.001…0.01]	0.0100	0.0100
ϖopt	[]	[1.5…10]	2.45	3.47	[1.5…10]	1.82	2.03
m	[1/s]	[5…20]	12.8	13.7	[5…20]	8.17	8.22
ν	[]	[2…8]	4.99	4.85	[2…8]	3.76	4.77
k^_hyd	[1/s]	[1…3]	NaN	1.45	[1…3]	NaN	2.27

Parameter values are given to three significant digits. The slack length of the serial elastic element constitutes for an exception with four digits due to the high sensitivity of the model with respect to this value (cf. .

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 () 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 L1 (absolute) and L2 (squared) distances for the purpose of showing consistency throughout the model variants.

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

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; green PCR) as well as the relative chemical potential of ATP (black: ). 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 (Fmax) as calculated by the optimization algorithm. The explanation of an increasing Fmax 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 Brel, 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 ΔUSEE, nll and ΔUSEE, 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 Fmax 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 Fmax. 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).

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−1 was found by an optimization routine, all other parameters were taken from Table 2. The unmodified and scaled force-length relation [, gray line] is given for reference. In (C,D) the procedure is mirrored for PLA, where κPLA = 0.165 s−1.

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 (DSDE = 10) and completely force independent (RSDE = 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 (Se = 4.5) in order to compensate a slow force decay at the end of the stimulation.

Maximum shortening velocity (Brel, 0/Arel, 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 Arel, 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 () seems to metabolize ATP, and thus accumulate phosphate, more rapidly than GAS (), 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 (L2) 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 Fmax 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 (L1) 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

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 104.

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; 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 () can be written as

which at full troponin-activity () results in a force decay rate of

see also Equations (37)–(39) for the latter proportionality. Indeed, the connection 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.

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 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 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).

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 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] with the characteristic time constants for condensation kcon and hydrolysis khyd. This ansatz together with the introduction of the state leads to the observed force decay (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 khyd 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 determines the corresponding characteristic time of decay. It is remarkable that 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/ 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 ADP2− and (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 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] 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] (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 khyd and kcon. 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 [Ca2+] sensitivity (Allen et al., 2008). The first component is associated with elevated [P] due to the breakdown of PCr. When [P] had been elevated from below 1 to ~14 mM, a force decrease of about 50% was observed in [Ca2+]-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] 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] 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 [Ca2+] sensitivity (Debold, 2016). It was shown that the force decayed by 10% while the myoplasmic [Ca2+] increased (Westerblad and Allen, 1991). This increase is interpreted as the result of reduced myoplasmic [Ca2+] buffering (Westerblad and Allen, 1993). Alterations of the myofibrillar [Ca2+] sensitivity might occur due to increase of elevated metabolites as [H+] and [P]. Nelson and Fitts (2014) observed at a pH of 6.2 in skinned muscle fibers (at 30°C) a decreased myofibrillar [Ca2+] sensitivity. The myofibrillar [Ca2+] sensitivity decreased even more when they added 30 mM [P]. It was concluded that both low pH and elevated [P] have a substantial effect on myofibrillar [Ca2+] 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] 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], which influences both actomyosin cross-bridge force generation and myofibrillar [Ca2+] 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 vmax in isometric contractions used for parameter fitting. Thus, uncertainty in parameter estimation is higher for vmax and Pmax compared to classic methods, in which vmax 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 (ΔUSEE, l ≈ 8 vs. 4.9% for GAS and 3.6% for PLA with classic methods), with a transition to the linear region at around Fmax and a less stiff behavior in the linear region (K ≈ 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 ≈ 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% Fmax; PLA: 9% Fmax) 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

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 Fmax ≈ 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%Fmax, which is in perfect agreement with the 15%Fmax 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

(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] 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).

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.

Abbreviation	Meaning
ADP	Adenosine diphosphate
AMP	Adenosine monophosphate
ATP	Adenosine triphosphate
Ca (Ca2+)	Calcium (ion)
CE	Contractile element
Cr	Creatine
GAS	M. gastrocnemius
MTU	Muscle tendon unit
ODE	Ordinary differential equation
Pi	Inorganic phosphate
PCr	Creatine phosphate
PEE	Parallel elastic element
PLA	M. plantaris
pH	Power/potential of hydrogen
SDE	Serial damping element
SEE	Serial elastic element

Symbol	Value	Unit	Meaning
ã	q~·μ~_ATP	[]	Activity
Arel, 0	Parameter	[]	Hill parameter
Brel, 0	Parameter	[1/s]	Hill parameter
DSDE	Parameter	[]	Shape parameter of the SDE force-force curve
ΔFSEE, 0	Parameter	[N]	Shape parameter of the SEE force-length curve
ΔUSEE, nll	Parameter	[]	Shape parameter of the SEE force-length curve
ΔUSEE, l	Parameter	[]	Shape parameter of the SEE force-length curve
ΔWasc	Parameter	[]	Shape parameter of the isometric force-length relation
ΔWdes	Parameter	[]	Shape parameter of the isometric force-length relation
F, FX	Free	[N]	Force, force of element X
Fe	Parameter	[]	Shape parameter of the eccentric force-velocity branch
F~_isom(ℓ_CE)	Length-dependent	[]	Normalized isometric force-length relation
Fmax	Parameter	[N]	Maximum isometric force at ℓCE, opt
F_PEE	Parameter	[]	Shape parameter of the PEE force-length curve
ΔG°	Reaction-dependent	[kJ/mol]	Standard-state free (Gibbs-)energy of reaction
K	Free or state-dependent	[]	Equilibrium constant
k	Free or state-dependent	[1/s]	(Chemical) reaction rate constants
k^_hyd	Parameter	[1/s]	ATP hydrolysis rate constant
ℓCE	State	[m]	Length of the contractile element
ℓCE, opt	Parameter	[m]	Optimal length of the contractile element
ℓ~_CE	ℓCE/ℓCE, opt	[]	Relative length of the contractile element
L_PEE	Parameter	[]	Shape parameter of the PEE force-length curve
ℓSEE, 0	Parameter	[m]	Slack length of the serial elastic element
M, mM, μM	100, 10−3, 10−6	[mol/l]	Molar, millimolar, micromolar
m	Parameter	[1/s]	Time constant (Hatze)
μATP, μATP, max	State-dependent	[kJ/mol]	(maximum) chemical potential of ATP
μ~_ATP	μATP/μATP, max	[]	Relative chemical potential of ATP
ν	Parameter	[]	Exponent (Hatze)
νasc	Parameter	[]	Shape parameter of the isometric force-length relation
νdes	Parameter	[]	Shape parameter of the isometric force-length relation
νPEE	Parameter	[]	Shape parameter of the PEE force-length curve
q~	q~_min≤q~≤1	[]	Troponin-activity (Hatze)
q~_min	≤ 0.01	[]	Minimum troponin-activity
RSDE	Parameter	[]	Shape parameter of the SDE force-force curve
Se	Parameter	[]	Shape parameter of the eccentric force-velocity branch
ϖ(ℓ~_CE)	ϖ_c·ℓ~_CE	[]	Length-dependent volumetric troponin C density
ϖopt	Parameter	[]	Constant (Hatze)
u	0 ≤ u(t) ≤ 1	[]	Neural stimulation
[X]	Free	[mol/l]	Molar concentration of substance X
[X]0	Free	[mol/l]	Initial molar concentration of substance X
[X]min	Free	[mol/l]	Minimum concentration of X