On the rate-limiting dynamics of force development in muscle Open Access

Skeletal muscle produces force many times slower than its neural excitation. The time course is due to activation dynamics that transform neural excitation into biochemical force-generating capacity, followed by the actual production of force – biophysically attributed to cycling of muscle cross-bridges (Fig. 1, center) – against series elastic tissue and the body skeleton. These dynamics are complicated by feedback, in that the resulting body motion also affects the force. The complexity of these interactions obscures which of the preceding processes are rate limiting and thus most responsible for the relative slowness of force production. This question could be addressed with quantitative muscle models, but there is little agreement or definition for the underlying dynamics. There are presently two main types of models for quantifying force production: (1) mechanistic cross-bridge models that focus on the biophysics of force production but not its influence on movements (Huxley, 1957; Zahalak, 1981; Campbell et al., 2018; Newhard et al., 2019; Liu et al., 2024), and (2) phenomenological ‘Hill-type’ muscle models (Hatze, 1977; Hof and Van den Berg, 1981a; Zajac, 1989; van Soest et al., 1993; Thelen, 2003; Mayfield et al., 2023) that interact with tendon and body but lack mechanistic dynamics. There is a need for integration and resolution of both biophysical and phenomenological approaches to mechanistically determine what limits skeletal muscle force development during movement.

Muscle force development is characterized by several types of experimental observations (Fig. 1). Whereas steady force production is described by well-known force–length (Gordon et al., 1966) and force–velocity (Hill, 1938) relationships (Fig. 1, ‘Steady response’, left), the dynamic effects are evident in time (Fig. 1, ‘Temporal response’, right) and frequency domains (Fig. 1, ‘Frequency response’, bottom). The temporal response shows how excitatory action potentials trigger relatively fast calcium release and binding, followed by relatively slow force development. The frequency response describes how force magnitude and timing vary as a function of excitatory input frequency (Partridge, 1965). For a muscle fiber, the action potential firing rate could be treated as input (Bigland-Ritchie et al., 1983; De Luca and Hostage, 2010) and its force output treated as output, characterized by a force-firing rate curve (Binder-Macleod and McDermond, 1992). For a whole muscle, the periodic electromyographic (EMG) signal could be treated as input (Hof and Van den Berg, 1981b; Zajac, 1989), with the output described by force magnitude (gain) versus the frequency of cyclic muscle contractions (referred to as contraction frequency). The force magnitude (gain) decreases with contraction frequency (Coggshall and Bekey, 1970; Calvert and Chapman, 1977; Crosby, 1978; Genadry et al., 1988; Hof, 1997; Winter, 2009) and its phase lag increases, similar to a low-pass filter (see Fig. 1, ‘Low-pass filter effect’, bottom). Although these effects have been well characterized experimentally for decades, it remains unclear how they relate to the underlying biophysical dynamics of muscle.

Force development is described by several types of muscle dynamics (Fig. 2). Musculoskeletal models (e.g. Zajac, 1989) usually distinguish between displacements of contractile and lumped elastic elements (Fig. 2A), between muscle dynamics for activation and contraction, and between muscle–tendon and body dynamics (Fig. 2B). Muscle dynamics are usually simulated with conventional Hill-type or cross-bridge muscle models, which include first-order activation dynamics and contraction dynamics. Activation dynamics transform neural excitation u(t) into activation a(t). Hill-type model contraction dynamics transform activation a(t) into muscle force F(t), phenomenologically mediated by the force–length and force–velocity dependencies (represented as feedback of muscle–tendon actuator length L_MT and velocity ⁠, Fig. 2B,C). For whole-body movements, it is usually necessary to simultaneously predict both muscle lengths and forces (e.g. L_M and F in Fig. 2C), which depend on dynamical interaction between body and series elasticity. Hill-type models typically accomplish this by treating muscle length as an internal state (Fig. 2B) and using elasticity to relate force and length changes (Zajac, 1989). This approach contrasts with more mechanistic models that focus on the dynamics of cross-bridge cycling (Fig. 1, center), which may be twice as slow as calcium activation (Zahalak and Motabarzadeh, 1997), and thus at least as important for force development. But most such models cannot be integrated with body and elasticity, except at impractically high computational cost (Fig. 2D; Lemaire et al., 2016; van Soest et al., 2019). Hill-type and cross-bridge muscle models have different states, dynamics and parameters, and cannot generally be interchanged within musculoskeletal models. It is also unclear whether either model, or perhaps a combination of the two, can explain force development for whole muscle.

A long-standing inconsistency is how models define muscle activation, sometimes referred to as active state. Hill (1949) originally treated active state as an inferred force capacity, not directly measurable but estimated from dynamic experiments (e.g. quick-release: Edman, 1970; Edman and Kiessling, 1971; Inbar and Adam, 1976). Others have provisionally defined activation mechanistically in terms of cellular calcium concentrations (e.g. bound to troponin, CaTrop: Ebashi and Endo, 1968; Hatze, 1977, 1978; Zahalak and Ma, 1990) but lacking direct measurements for such concentrations. Most models have therefore defined activation more loosely, as arising from first-order activation dynamics that, when scaled by the steady force–length and force–velocity relationships (Fig. 1, left) and interfaced with contraction dynamics, yield the muscle force output (Zajac, 1989; Sandercock and Heckman, 1997; Perreault et al., 2003). Associated activation (and deactivation) time constants are determined ad hoc in some models (Zajac, 1989; Thelen, 2003; Millard et al., 2013), and tuned to match overall force development in others (Curtin et al., 1998; van Zandwijk et al., 1998; Lemaire et al., 2016). Consequently, estimated time constants differ considerably, up to 10-fold between definitions (Calvert and Chapman, 1977) and severalfold between similar-seeming Hill-type models (e.g. van Zandwijk et al., 1998; Millard et al., 2013). Another approach is to fit a muscle's EMG-to-force transfer function (Calvert and Chapman, 1977; Lloyd and Besier, 2003; Buchanan et al., 2004; Lee et al., 2011, 2013). A transfer function allows for more complex dynamics to non-parametrically describe force production across a range of excitation frequencies, but is also phenomenological and cannot indicate where the dynamics occur. Whereas the early mechanistic approaches to activation (e.g. Hatze, 1977, 1978) have lacked data, the more empirical approaches have lacked mechanisms.

However, greater mechanistic insight is enabled by more recent studies. Fast-acting calcium dyes have enabled higher temporal resolution of calcium binding dynamics (Zhao et al., 1996), which play a substantial role in energy expenditure (Barclay et al., 2007). Other studies have revealed additional processes in force production (Irving, 2017). Whereas calcium activation (e.g. CaTrop) was previously thought to initiate cross-bridge cycling, more recent studies have found significant intermediate dynamics, such as for movement of tropomyosin (Kress et al., 1986; Walcott, 2014) that dynamically affects cross-bridge binding site availability on the thin filament (McKillop and Geeves, 1993; Campbell, 2014), and for refolding of myosin heads out of a ‘super-relaxed’ state to enable subsequent cross-bridge binding (Reconditi et al., 2011; McNamara et al., 2014; Linari et al., 2015; Campbell et al., 2018; Nag and Trivedi, 2021; Liu et al., 2024). These processes can take significant time (e.g. 20–50 ms; Kress et al., 1986; Reconditi et al., 2011), perhaps exceeding that for calcium release and binding to troponin (Hollingworth et al., 1996; Baylor and Hollingworth, 2003). Such dynamics have primarily been studied at the cellular and molecular level thus far (Irving, 2017; Campbell et al., 2018), and could be rate limiting for whole muscles as well.

Here, the primary objective was to quantify the rate-limiting dynamics of muscle force development. Our approach begins (Part I) by examining prior experimental data for calcium activation and force development dynamics. The data superficially appear to follow first-order processes, but closer examination will show that the activation dynamics seem to vary based on the stimulus train. This suggests the need for quantitative modeling to test whether a more self-consistent set of dynamics can explain such data. We therefore (Part II) apply conventional Hill-type and cross-bridge muscle models previously employed in movement simulation to these data, and show that no such conventional model can self-consistently explain muscle force development well, because a new set of dynamics is needed. Finally (Part III), we propose an integrative muscle model to rectify this disparity. We integrate calcium activation with two types of contraction dynamics combining simplified cross-bridge dynamics with a Hill-type interface to tendon and body dynamics. To this we add a new set of ‘force facilitation’ dynamics, derived from Part II. The model is tested with a single set of self-consistent parameters, and will be shown to reproduce all three categories of empirical muscle responses (Fig. 1). The model shows that the proposed force facilitation dynamics are necessary to explain relatively slow muscle force development.

Biological muscle exhibits distinct patterns of calcium activation and force development (Fig. 3). Unavailable at the time of Hill (1949), simultaneous recordings of calcium transients and isometric force development now reveal that both resemble first-order processes, but with markedly different time constants. Calcium imaging studies (employing fluorescent dyes in vitro) show that calcium transport and binding to troponin are relatively fast (Baylor et al., 1983; Zhao et al., 1996), suggesting that other processes are necessary to explain slow force development. These data present an opportunity to test muscle models more quantitatively than previously possible.

Two types of responses are examined here, for step tetanic stimulus (Fig. 3A,B) and single twitch (Fig. 3C,D). Tetanic calcium activation (CaTrop from mouse; Hollingworth et al., 1996) is quite fast (Fig. 3A). Ignoring very fast ripples for individual action potentials that do not appreciably affect the force transients (Fig. 3B), the remaining overall time course resembles a saturating exponential, with a 90% rise time (τ₉₀, about 2.3 times the exponential time constant) of about 5 ms (Fig. 3A, first-order activation). The resulting force development transient also resembles a saturating exponential but more than 20 times slower, with a τ₉₀ of 115 ms (Fig. 3B, first-order force). Additional insights are revealed by the twitch response, which highlights the slower calcium deactivation transient, about 6 times slower than activation (Fig. 3C). The corresponding twitch force rises and decays much slower than calcium activation (Fig. 3D).

Although these observations all superficially resemble first-order dynamics, they are unlikely to be explained by Hill-type models. Isometric conditions are illustrative because they eliminate Hill-type contraction dynamics, thus isolating first-order activation dynamics as the sole explanation available for force development (e.g. Hatze, 1977; Zajac, 1989; Winters, 1995; Thelen, 2003). It thus appears that Hill-type models should be able to reproduce either fast calcium activation or slow force development, but not both simultaneously. In contrast, cross-bridge models have additional dynamics, not present in Hill-type models, that might help account for slower force development. We therefore next examine how both types of conventional muscle models (i.e. Hill-type and cross-bridge) compare against these data.

We quantitatively tested whether conventional Hill-type and cross-bridge muscle models (Fig. 2) could reproduce empirical activation and force development dynamics. We consulted two data sources: (1) in vitro mouse data for calcium activation and force development as above, and (2) in vivo human data of force development (see Supplementary Materials and Methods for details). The latter lack calcium activation estimates, but do include two types of force responses, and thus two ways to test the models. Because of the wide range of model definitions for activation, we considered two extremes for the first-order time constant: a faster and more mechanistic one consistent with calcium-based activation, and a slower and phenomenological one matching force data. Although the data are available in the literature, neither Hill-type nor cross-bridge models have previously been tested integratively for both activation and force development. Below, the model dynamics are briefly described and then tested against data.

We implemented both Hill-type and cross-bridge models, with some features in common. They share common lumped elements for parallel elasticity, series elasticity (attributed mainly to tendon) and contractile element (Fig. 2A, PE, T and CE, respectively; CE and PE share the same length L_M). Series elasticity (T in Fig. 2A) is modeled with an exponential-to-linear relationship (Winters, 1995), and parallel elasticity with a quadratic relationship (van Soest and Bobbert, 1993). The active force–length relationship is modeled with a quadratic dependency on length L_M (van Soest and Bobbert, 1993) to resemble an inverted U-shape (Bohm et al., 2018; Bohm et al., 2019) (see Table S1 for parameter values). Both Hill-type and cross-bridge models therefore share the same parameter values for these mechanical features.

We use first-order dynamics (Fig. 2B) between muscle excitation u and activation a, as previously used in both Hill-type and cross-bridge muscle models (Hatze, 1977; Hof and Van den Berg, 1981a; Zajac, 1989; van Soest and Bobbert, 1993; Sandercock and Heckman, 1997; Thelen, 2003; Perreault et al., 2003; Lemaire et al., 2016; Mayfield et al., 2023). We adopted separate time constants for activation and deactivation (i.e. faster forward and slower backward reactions) because of their different underlying mechanisms (i.e. calcium release and uptake, respectively, see Appendix for details). As with most existing muscle models, motor unit recruitment and firing rate are lumped into a single quantity u, summarizing excitation for the entire muscle.

The two extremes for muscle activation dynamics are as follows. We define ‘calcium-based activation’ physiologically as the relative amount of calcium bound to troponin regulatory sites (CaTrop) (Ebashi and Endo, 1968; Hatze, 1977; Zahalak and Ma, 1990; van Soest and Bobbert, 1993; Lemaire et al., 2016), based on empirical data (Hollingworth et al., 1996). In contrast, ‘force-based activation’ refers to a latent and not directly measurable state, selected to approximately agree with force development transients (Zajac, 1989; Sandercock and Heckman, 1997; Umberger et al., 2003; Thelen, 2003; Perreault et al., 2003). Only one model parameter was varied between definitions, i.e. the (forward) activation dynamics time constant. An implication is that a calcium-based activation will reproduce calcium concentration data but not necessarily forces, and a force-based activation the converse.

We next examine how these muscle models compare with empirical activation and isometric force. Mouse and human comparisons are presented in separate figures (Figs 4 and 5), each with both step tetanic excitation and twitch excitation (top and bottom rows, respectively), simulated with two extremes of activation time constants: fast calcium-based (Figs 4A and 5A) and slower force-based (Figs 4B and 5B) as subfigures. Within each subfigure, the comparisons include calcium activation (CaTrop), where available, and force development. As with previous modeling studies, twitch and tetanus excitations were simulated using a 5 ms unit pulse (Thelen, 2003; Mayfield et al., 2023) and unit step input, respectively (Curtin et al., 1998; Lemaire et al., 2016).

These comparisons show that activation dynamics contribute more to force production than other muscle dynamics. This is because (near-) isometric conditions downplay the effect of contraction dynamics and interaction with body, and cross-bridge cycling dynamics needed to be relatively fast to comply with force–velocity relationships (see Appendix for details). As the two models share the same activation dynamics, both Hill-type and cross-bridge models produce very similar simulations (dotted versus dashed blue traces, respectively in Figs 4 and 5). We therefore describe general ‘model’ results here, only distinguishing between Hill-type and cross-bridge models where specifically needed.

We first compared the models with in vitro mouse data from fast-twitch extensor digitorum longus muscle (Fig. 4), including simultaneous recordings of CaTrop and force (Hollingworth et al., 1996) (see Supplementary Materials and Methods for details). As expected, a calcium-based activation model can successfully replicate CaTrop transients for both step tetanic and twitch excitation (Fig. 4A). However, combining this activation with contraction dynamics results in forces far faster than empirical data. Additionally, the models' twitch force response exhibits far greater amplitude than data show (Fig. 4A). Alternatively, selecting a slower force-based activation time constant enables models to reproduce step tetanic force development (Fig. 4B). However, the associated activation dynamics must then be far too slow to match empirical CaTrop data. Moreover, a model proficient in replicating step tetanic force performs inadequately for twitch force, with far too low an amplitude (Fig. 4B). Consequently, conventional models can predict either activation or force transients observed experimentally, but not both (Fig. 4). Altogether, the contraction dynamics of both Hill-type and cross-bridge muscle models prove to be excessively fast to account for these empirical data.

Similar results were obtained for in vivo human quadriceps muscle data (Fig. 5), including dynamometer-based estimates of muscle force (Verges et al., 2009) (see Supplementary Materials and Methods for details). Lacking calcium activation data for humans, the models are tested by their ability to reproduce force development data for both step tetanic and twitch stimuli. Again, we consider two possibilities for model activation time constants, selected to correspond to (assumed) calcium-based activation or to agree with empirical force data. As expected, calcium-based activation is relatively fast, and the subsequent model contraction dynamics are insufficient to explain the relatively slow forces observed experimentally. The models produced far faster tetanic forces and larger twitch forces compared with empirical data (Fig. 5A). The alternative of slower, force-based activation was able to reproduce the step tetanic force response well, but then greatly underestimated twitch force amplitude (Fig. 5B). Altogether, even in the absence of empirical CaTrop activation data, current first-order models of activation dynamics cannot reproduce both twitch and tetanic force responses.

These findings confirm the suspicions from Part I, and show that conventional muscle models cannot simultaneously predict twitch and tetanic force development data. Models with calcium-based activation can reproduce fast calcium activation dynamics (when empirical data are available, Fig. 4), but not the observed slow force development, because their contraction dynamics are too fast to reproduce data. Models with force-based activation can reproduce tetanic force data (Figs 4 and 5), but at the expense of not agreeing with twitch force data (Fig. 5), and/or requiring unrealistically slow calcium activation (Fig. 4). This evidence suggests that muscle exhibits an additional set of dynamics, not explainable by any first-order model of activation dynamics or by either type of contraction dynamics.

We next propose a new muscle model that resolves many of the above issues, called the CaFaXC (calcium-facilitation–cross-bridge-contraction) model (Fig. 6). It integrates several features of existing models, including calcium-based activation dynamics (c, representing physiological CaTrop; Fig. 6A) and contraction dynamics that encompass both length and body feedback and cross-bridge cycling dynamics. However, it also inserts a second and slower set of first-order dynamics, termed force facilitation dynamics (Fig. 6B), before contraction dynamics. Studies have revealed force development is affected by processes such as tropomyosin movement (e.g. Moore et al., 2016) and myosin unfolding (e.g. Linari et al., 2015), although it is presently uncertain which of them is rate limiting. Rather than modeling them explicitly, we provisionally model their lumped effect with an additional latent force facilitation state r that reflects relatively slow force development. Below, we first describe the model and then test it against experimental calcium and force data (Fig. 1).

The model's parameters were determined as follows (see Appendix for details). The calcium-based activation dynamics have activation and deactivation time constants determined either from empirical CaTrop (for comparison with mouse data, described below) or from twitch and tetanic force responses (for human, described below). Force facilitation also has forward and backward time constants, here tuned to match empirical force, similar to force-based activation in Hill-type models. The contraction dynamics require four cross-bridge rate parameters, here constrained by the force–velocity relationship, similar to that described previously (Huxley, 1957; Zahalak, 1981; van Soest et al., 2019). Three of them may be mapped to three Hill-type force–velocity parameters; for example, maximal shortening velocity, force–velocity curvature and eccentric force asymptote (see Appendix). The fourth parametric degree of freedom is a free parameter, unexamined here, that effectively relates shortening velocity to energetics, somewhat like some models that link Hill-type models to phenomenological energetics estimates (e.g. Umberger et al., 2003). The proposed model also has length and force scaling parameters to ensure compatibility between cross-bridge strains and muscle length changes (Lemaire et al., 2016). These scaling parameters may be mapped from optimal fiber length and maximum isometric force of Hill-type models (see Appendix). The model's parameters may therefore largely be mapped from Hill-type models, except for two additional time constants for force facilitation dynamics.

The CaFaXC model was tested against step tetanic and twitch temporal response data (Fig. 7). Starting with mouse tetanic response (Fig. 7A), the fast activation dynamics yield fast activation rise (τ₉₀=5 ms), in line with empirical CaTrop data (Hollingworth et al., 1996). Using the empirical force development data to determine the force facilitation dynamics, the model also reproduces relatively slow force development (τ₉₀=102 ms, Fig. 7A). These time constants provide quantitative contributions to overall force development time: calcium activation 5%, force facilitation 76%, cross-bridge dynamics 3% and contraction dynamics 16%. Without modification, the same model reproduces a twitch response that shows fast activation and slow force. Unlike conventional models (Fig. 4B), the integrated model predicts much smaller twitch force amplitudes compared with tetanic force (1:5 ratio, Fig. 7A), similar to experimental data. Both tetanic and twitch responses were dominated by force facilitation dynamics and reproduced data well.

The model also compares well with human data (Fig. 7B). Lacking human CaTrop measurements comparable to those of mouse, the model was instead compared against the combination of twitch and tetanus data. With a single set of parameters, the model reproduces both of the relatively slow twitch and tetanic force responses (τ₉₀=199 ms, Fig. 7B) reasonably well. The contributions were calcium activation 3%, force facilitation 67%, cross-bridge dynamics 3% and contraction dynamics 27%. Again, it produces a small twitch-to-tetanus ratio similar to data (1:5, Fig. 7B), unlike conventional models (Fig. 5B). Although the associated activation dynamics are not presently testable for humans, the model predicts relatively fast activation (τ₉₀=5 ms, Fig. 7B) and slow force facilitation, consistent with mouse data. Apart from different time scales, the model suggests that activation and force development dynamics are qualitatively similar for mouse and human. A self-consistent model for each can simultaneously reproduce activation (in mouse) and force responses (in mouse and human), including both twitch and tetanus.

We next examined whether a single set of CaFaXC parameters for human quadriceps muscle can reproduce multiple empirical observations. First, as the model's steady responses were constrained to empirical force–length and force–velocity relationships for human quadriceps, these are reproduced relatively well (R²=0.85 and R²=0.96, respectively; Fig. 8, left). The model also reproduces both twitch and tetanus force development (R²=0.95 and R²=0.98, respectively; Fig. 8, top middle and right), with tetanus rise time and twitch amplitude similar to empirical values. The same model then makes independent predictions of frequency response data, including the EMG-to-force relationship as a function of cyclic contraction frequency (Fig. 8, bottom right). Using the same parameter values as above, the model shows a decrease in gain and an increase in phase lag between EMG and force, in agreement with empirical data (R²=0.93 and R²=0.82, respectively). In addition, the model predicts a saturating increase for the force–firing rate relationship, similar to experimental data (R²=0.96; Fig. 8, bottom middle). Altogether, the model self-consistently reproduces multiple muscle responses, including classic steady force–length and velocity relationships, temporal responses and frequency responses, across a range of conditions.

This study aimed to quantitatively determine which dynamics underly muscle force development. Relatively slow dynamics are apparent from muscle's temporal and frequency responses, but are not captured by conventional muscle models. We found that an additional process, termed force facilitation, is needed to reconcile slow force development with fast calcium activation. We next examine the relative timing contributions from our integrated model, consider the mechanisms underlying the proposed force facilitation dynamics and discuss the biological relevance and implications for movements.

These results resolve long-standing questions regarding activation dynamics. We found that any first-order definition of activation, whether coupled to Hill-type or to cross-bridge mechanics, can only reproduce one of either tetanic activation or force development (Fig. 4), or one of tetanic and twitch forces (Figs 4 and 5). With calcium imaging (Baylor and Hollingworth, 2003; Hollingworth et al., 1996; Zhao et al., 1996), activation can be defined specifically, for example as CaTrop here, or as free or total calcium concentration. But to agree with such data, a calcium-based definition requires fast activation dynamics. We found a time constant (2 ms) to be compatible with both human and mouse data, in agreement with fast calcium dynamics previously reported across a range of species and muscles (Baylor and Hollingworth, 2003; Rome et al., 1996). Such fast calcium dynamics therefore can account for only 3–5% of overall force development time. This is an order of magnitude faster than in previous models that have qualitatively attributed activation dynamics to calcium processes (e.g. van Soest et al., 1993; Thelen, 2003; Mayfield et al., 2023). Both data and model now show that calcium dynamics, while critical to muscle function, are too fast to account for actual force development.

Additional dynamics are necessary to explain slow force development in muscles. In most biomechanics studies, muscle dynamics have comprised some subset of calcium-based activation, cross-bridge mechanics and Hill-type contraction against series elasticity and body. These dynamics can be resolved by applying empirical data and enforcing dynamical constraints. For example, Huxley's (1957) cross-bridge mechanics both explain and are constrained by the muscle force–velocity curve. This constraint required quite fast rate constants here, so that cross-bridge dynamics accounted for only about 3% of force development time. Hill-type contraction dynamics also contribute to force development as a result of the force–velocity curve, as muscle fibers shorten against series elasticity (despite isometric conditions for the limb). Constrained by experimental data for human quadriceps, such contraction dynamics accounted for a substantial 27% of overall force development time. However, even the combination of calcium activation with conventional Hill-type and cross-bridge contraction dynamics can only explain about one-third of overall force development time.

We have therefore introduced a new latent state, provisionally termed force facilitation, to account for the missing force development time. Specifically, a first-order process following calcium dynamics is sufficient for compatibility with both calcium and force data. It presently lumps together more mechanistic processes (Fig. 6C) with a single time constant empirically determined to be quite slow, accounting for most (about 76% for mouse and 67% for human) of force development time. The need for such dynamics may not have been appreciated previously, because few models have addressed calcium activation and force development directly, and have thus not exposed the temporal gap between the two. The proposed force facilitation is necessary to account for that gap, and appears to be the time-dominating process in muscle's ability to produce force, also known as latent activation.

The integrated model also provides insight regarding the frequency response behavior of whole muscle (Figs 1 and 8). Muscle has long been observed to resemble a low-pass filter of excitation (Partridge, 1965; Stein et al., 1972; Milner-Brown et al., 1973; Bawa and Stein, 1976; Crosby, 1978; Winter, 2009). For example, force magnitude per unit excitation (gain) decreases with contraction frequency, similar to a critically damped second-order system with a relatively slow natural frequency (e.g. 2–3 Hz; Bawa and Stein, 1976; Milner-Brown et al., 1973). Muscle is thus far slower than explained by calcium activation (Fig. 7) or latent activation dynamics of typical Hill-type models (e.g. Zajac, 1989; Thelen, 2003; Millard et al., 2013), but the discrepancy was unclear because frequency response methods and EMG-to-force transfer functions (e.g. Calvert and Chapman, 1977; Lloyd and Besier, 2003; Buchanan et al., 2004; Lee et al., 2011, 2013) usually provide little physical mechanism or location for the identified time constants. Our model with more mechanistic calcium and contraction dynamics suggests that force facilitation may be the dominant (slowest) dynamics underlying the low-pass filter effect (Fig. 7), followed by contraction dynamics.

Force facilitation dynamics may be relevant to modeling of human movement. Lacking rate-limiting force development dynamics, Hill-type models with ad hoc first-order activation are often too fast. In musculoskeletal simulation (e.g. OpenSim: Seth et al., 2011), such models have been observed to underestimate the amplitude of muscle excitation transients for faster walking speeds (Luis et al., 2022), and EMG-to-force delays for walking (Luis et al., 2022) and hopping (Jessup et al., 2023). Such models also underestimate the energetic cost of faster movements; for example, only predicting 7% and 40% of the cost increases for cyclic torque production (van der Zee and Kuo, 2021) and reaching movements (Wong et al., 2021). Other models have used slower activation time constants that agree better with (tetanic) force development data (Curtin et al., 1998; Lemaire et al., 2016). But such models may underestimate the amplitude of twitch forces (see Figs 4 and 5), or fail to represent higher order dynamics evident in reaching movements (Murtola and Richards, 2023). In addition, attribution of timing to calcium rather than force facilitation dynamics could lead to incorrect estimates of calcium transport and its energetic cost (Szentesi et al., 2001; Barclay et al., 2007), potentially leading to underestimated energetics as stated above. The inclusion of force facilitation dynamics, whether in the present model or otherwise-conventional Hill-type models, may improve how well musculoskeletal simulations can perform when tested against independent experimental data. Such models can potentially improve upon predictions of muscle excitation, forces and energetics during movement.

The present model also integrates together previously disparate dynamics. To date, Hill-type models have generally included dynamical interactions with body and series elasticity but not cross-bridge dynamics, and cross-bridge models the converse. Their combination has previously been found to be computationally impractical (van Soest et al., 2019), which we addressed by adopting the distribution moment approach (Zahalak, 1981), resulting in modestly greater complexity than Hill-type models (see Appendix). We found both types of dynamics to be relatively fast but nonetheless important for modeling movement. Cross-bridges are the primary explanation for the force–velocity relationship, which must arise from the same mechanics that govern force production. Such mechanics also predict transient muscle phenomena such as short-range stiffness (Campbell and Lakie, 1998), observable in upright balance (De Groote et al., 2017) and implicated in proprioception (Simha and Ting, 2024). Cross-bridge cycling also mechanistically predicts the energetics of contraction (Ma and Zahalak, 1991). This contrasts with Hill-type models that, because of phenomenological curve fits, enforce no mechanistic consistency (e.g. between force production and force–velocity) and have no mechanism for short-range stiffness. Such models are thus susceptible to overfitting (e.g. Figs 4 and 5), making it difficult to make generalizing predictions. Here, we guarded against overfitting by using twitch, tetanus and/or calcium data to identify force facilitation, and then independently testing both the low-pass filter effect and the force–firing rate relationship (Fig. 8). The proposed model integrates cross-bridge dynamics for mechanistic consistency into a dynamical interface with the body derived from the Hill-type approach.

Although we advocate for force facilitation dynamics, its underlying mechanisms remain unclear. Potential processes include cooperative interaction between thin and thick filaments (Campbell, 2014; Walcott, 2014) and facilitative recruitment of myosin heads (Reconditi et al., 2011; McNamara et al., 2014; Linari et al., 2015; Irving, 2017; Nag and Trivedi, 2021; Ma et al., 2023), mediated by other protein interactions (Song et al., 2021; Squarci et al., 2023; Hessel et al., 2024). We provisionally lumped an overall effect into the force facilitation state r, identifying its time constant from force development data. This latent state is presently the least mechanistic feature of our model, and might seem like little improvement over Hill's (1949) original active state, but force facilitation is more easily identifiable from experiments (e.g. Fig. 7), and more clearly localized to the processes between calcium activation and actual cross-bridge binding. It suggests that Hill's active state might be mechanistically interpretable as a measure of cross-bridge binding, such as a fraction of attached cross-bridges. Further experimentation and modeling is needed to resolve how muscle force is actually facilitated (e.g. Longyear et al., 2017; Campbell et al., 2018; Liu et al., 2024; Squarci and Campbell, 2024). There is a need to assess and compare force facilitation across slow to fast muscle fiber types, between species, and within a species but between muscles. But in contrast to the Hill-type models that also reproduce tetanic force development data (Curtin et al., 1998), force facilitation here highlights the timing gap between calcium activation and force development.

Additional limitations of the CaFaXC model include the need for two additional model parameters compared with Hill-type models. The model's contraction dynamics parameters may be mapped directly from Hill-type models (see Appendix), but the separate calcium and force facilitation dynamics require two additional parameters, for forward and backward calcium-based activation rates, distinct from current Hill-type models. These parameters are identifiable from in vitro CaTrop estimates (Fig. 7A) where available. But for human, we used (in vivo electrically stimulated) force responses (Fig. 7B) to provisionally identify both activation and force facilitation parameters simultaneously. This experimental procedure is in principle applicable to a variety of human muscles, facilitating quantitative determination of force development that contrasts with the ad hoc parameters used in many Hill-type models. However, a disadvantage is that the model's calcium activation parameters have yet to be tested against independent measures of CaTrop for human, presently unavailable. Such data would also be helpful for testing more complex models of calcium transport and binding (e.g. Baylor and Hollingworth, 2003; Kim et al., 2015). We consider the CaFaXC model to represent the minimum necessary to reproduce both calcium activation and force development data while having identifiable parameters.

We have thus far used relatively simple cross-bridge dynamics, as in Huxley's classic model (Huxley, 1957). While more complex cross-bridge models have become available (Campbell and Lakie, 1998; Campbell and Moss, 2002; Campbell, 2014; Walcott, 2014; Campbell et al., 2018; Kosta et al., 2022; Liu et al., 2024), the classic model has remained valuable (Lemaire et al., 2016; Newhard et al., 2019; van Soest et al., 2019), in part because it mechanistically explains force–velocity and energy–velocity relationships (Huxley, 1957). These relationships serve as constraints on the model, which must also explain isometric force and force development transients and temporal and frequency response dynamics (Fig. 8). This requires only a small number of parameters and helps guard against model over-fitting. While other, more complex, cross-bridge models are valuable for detailing features such as weakly bound cross-bridge states (Huxley and Simmons, 1971; Jarvis et al., 2021), and force transients during lengthening (Lombardi and Piazzesi, 1990), we found the classic cross-bridge model sufficient to test force development dynamics, and consider it adaptable to more complex models.

The present model is amenable to additional features such as force transmission and energetics. Musculoskeletal models often include muscle pennation as part of force transmission to the skeleton (Thelen, 2003; Millard et al., 2013; D'Hondt et al., 2024), and some also estimate metabolic energy expenditure (e.g. Umberger et al., 2003), albeit phenomenologically. The present model's contractile apparatus has the same external interface as Hill-type models, and is thus readily adaptable to pennation. Owing to the DM approximation, CaFaXC's modest computational cost (about twice that of the Hill-type model, see Appendix) facilitates implementation in movement simulation. The model is also adaptable to energetics, because it includes cross-bridge cycling and calcium activation, which are mechanistically implicated in energetics (Szentesi et al., 2001; Barclay et al., 2007). It remains to be tested whether the CaFaXC model predicts the energetics better than Hill-type models, but we consider the mechanistic approach to have greater potential, particularly because most muscle experiments are aimed at mechanistic understanding rather than phenomenological curve fitting.

The rate-limiting dynamics of force development are due to neither calcium activation nor cross-bridge cycling, nor to shortening against elastic tissues, which are all too fast to explain experimental data. Rather, an intermediate set of ‘force facilitation’ dynamics appear necessary to explain time- and frequency-dependent aspects of muscle force development in a self-consistent manner. A model of force facilitation introduces a relatively slow, dynamic relationship between muscle excitation and force, with potential applications to predicting the mechanics and energetics of movements.

We modeled the human quadriceps as a single muscle–tendon, lumping properties of the four individual heads together. Parameter values were either obtained from previous models or derived from previously reported data (see Table S1). Parameter values for mouse were either scaled from human quadriceps or derived from previously reported data (see 'Parameter values', ‘Parameter values for mouse extensor digitorum longus muscle', below).

The contractile element force–length relationship was assumed to be a quadratic function of muscle fiber length f_L(L_M) (van Soest and Bobbert, 1993), specified by maximal isometric force (F_max), muscle fiber optimum length (L_M,opt) and a width parameter, respectively (see Table S1). For the Hill-type model, we employed a standard force–velocity relationship (Umberger et al., 2003) specified by maximal shortening velocity v_max, curvature and eccentric force asymptote (see Table S1). Parameter values were chosen to approximately agree with torque–angular velocity data for human knee extension (Westing et al., 1990). For the cross-bridge model and CaFaXC model, we fitted attachment and detachment rate constants to yield the Hill-type force–velocity relationship in steady-state conditions (see ‘Parameter estimation’, ‘Cross-bridge cycling rates’, below).

The series elastic element force was assumed to increase linearly with strain for large forces, with an exponential toe region at small forces as described previously (Winters, 1995). Given an estimated tendon slack length (van Soest et al., 1993), we chose the remaining parameters (see Table S1) to agree with previously reported ultrasonography-based series elastic length changes in human quadriceps (van der Zee and Kuo, 2021).

The parallel elastic element force was assumed to increase quadratically with strain above passive element slack length, and to reach maximal isometric force F_max at a specified strain as described previously (van Soest et al., 1993), with parameter values from the literature (Thelen, 2003).

When incorporating cross-bridge dynamics into standard muscle modeling geometry (Fig. 2A), one needs to ensure that the forces in the three elements (i.e. CE, PE, T) are compatible. This is trivial in Hill-type models, because CE force F_CE is determined from forces in the elastic elements (i.e. F_PE and F). This is not trivial in cross-bridge models, because F_CE must simultaneously be compatible with forces in elastic elements and with contractile element cross-bridge dynamics. As conceived previously (Stoecker et al., 2009; Lemaire et al., 2016), elastic compatibility is accomplished through applying a constraint on the time-derivatives of the forces (i.e. ⁠) at each integration step.

As described in the main text, the activation dynamics time constants of the conventional muscle models were fitted to approximately agree with either empirical calcium activation or force development. The proposed CaFaXC model has one additional set of time constants, for a total of four time constants. These time constants were chosen to best match four features of force development: (1) tetanus rise time (τ_act, τ_fac), (2) twitch (or tetanus) decay time (τ_deact, τ_defac), (3) twitch time to peak (τ_act, τ_deact, τ_fac) and (4) twitch amplitude (τ_act, τ_deact, τ_fac). The additional experimental requirement for parameter estimation of the proposed model compared with conventional models is therefore assessment of twitch forces. Fortunately, twitch and tetanus require very similar experimental methods, and typically only differ in the duration of stimulation (brief or long). Here, we chose these time constants to approximately agree with twitch and tetanic force development for human quadriceps, while also evaluating force development in other conditions (e.g. force–firing rate, force–contraction frequency). The values for the activation dynamics time constants are in line with previous estimates of active state rise and decay in human muscle (Hof and Van den Berg, 1981b). Time constants are reported in the main text (see figure legends).

Parameter values for mouse extensor digitorum longus (EDL) fast-twitch fibers were estimated for the in vitro conditions considered here (Hollingworth et al., 1996). Hollingworth and colleagues (1996) estimated that 1–10 fibers were producing force, and that the average fiber diameter was 43 µm. Assuming an average of 6 isolated fibers, and a maximal tetanic tension of 390 kN m⁻² (Lännergren and Westerblad, 1987), maximal force of this fiber bundle was estimated at 0.0017 N. Muscle optimum length L_M,opt was estimated at 0.014 m (Hakim et al., 2013) and all other length parameters were assumed to scale linearly with the human quadriceps values. Maximal shortening velocity v_max was estimated at 10 L_M,opt s⁻¹, considering that: (1) mouse EDL is believed to be uniformly fast twitch (Hollingworth et al., 1996) and (2) smaller mammals typically have faster maximal shortening velocities than larger mammals (Wakeling et al., 2012). Because activation-dependent v_max scaling is thought to arise from recruitment of faster motor units (Petrofsky and Phillips, 1981), we did not apply this to the fast-twitch mouse muscle fibers (i.e. μ=1). Force–velocity curvature was estimated at 0.3, deemed appropriate for fast-twitch fibers (Wakeling et al., 2012). As with the human quadriceps, mouse EDL cross-bridge cycling rates were chosen as to fit the force–velocity relationship. The same value of h was used for cross-bridge scaling, but a smaller sarcomere length of 2.4×10⁻⁶ m (Goulding et al., 1997) yielded a slightly larger geometric scaling parameter γ (about 0.010).