# Muscle Physiology and Modeling

 George A. Tsianos and Gerald E. Loeb (2013), Scholarpedia, 8(10):12388. doi:10.4249/scholarpedia.12388 revision #143500 [link to/cite this article]
Post-publication activity

Curator: Gerald E. Loeb

Because they are so important for survival, muscles have evolved many complex features that enhance their strength, speed, efficiency and controllability. The model presented here consists of quantitative functions that represent each of the anatomical structures and physiological processes that underlie such features in typical mammalian skeletal muscles. These include the processes of recruitment, activation, active and passive force generation and energy consumption for the full physiological range of kinematics for both slow- and fast-twitch motor units. By selecting appropriate parameters, the model can be made to represent any specific normal or pathological muscle.

## Introduction

Quantitative models of muscle contraction are crucial for understanding neural control of movement. The nervous system must coordinate commands to a large set of muscles in a precise sequence. Characterizing their force producing capabilities provides insight into the feasible set of muscle activation patterns that lead to a desired behavior. The actual strategy selected from this set will likely depend on other factors such as metabolic energy consumption, whose minimization is clearly advantageous for survival. Both force generation and energy expenditure depend complexly on the commands from the nervous system, muscle length, and the rate at which muscle length changes. Computational models that relate these relationships to the underlying anatomy provide a means to validate theories of muscle physiology. Computational models that accurately capture these relationships are an essential component of complete models of musculoskeletal systems (see Musculoskeletal Mechanics and Modeling).

Neural control of muscles arises from a combination of commands from the brain as well as feedback from the periphery that are integrated by interneuronal circuits in the spinal cord. Muscle itself is a major source of the peripheral feedback, as it possesses many specialized mechanoreceptors that are sensitive to stretch and tension (see Proprioceptors and Models of Transduction). Both the stretch and tension experienced by muscle depend on the motion of the skeletal segments to which they attach; for a detailed explanation of the interactions among muscles, the skeleton, and environmental objects, see Musculoskeletal Mechanics and Modeling.

The physiological properties of muscle can change over time. In the short term, force generation may increase or decrease as a result of chemical changes associated with potentiation and fatigue, respectively. In the longer term, muscle morphology and physiology may change as a result of trophic responses to patterns of use. All of these effects tend to be specific to the various muscle fiber types, making it important to develop muscle models that reflect the subpopulations of fiber types and the relative recruitment of the motoneurons that control them.

## Relevant physiology

Each muscle is controlled by a group of motoneurons known as a motor pool or motor nucleus. All motoneurons in such a pool generally receive the same drive signals, although there are some exceptions (Loeb and Richmond, 1989). Each motoneuron along with all of its muscle fibers is defined as a motor unit, whose recruitment and firing rate in response to the same drive varies substantially depending on the size and impedance of the motoneuron (see Frequency-recruitment). Larger motoneurons tend to innervate larger numbers of larger diameter muscle fibers.

The force production and corresponding energy consumption of each motor unit depend mainly on its firing rate, muscle fiber length and velocity. Because muscle fibers within a motor unit generally have the same contractile properties, they can be lumped into one mathematical entity whose force generating capacity is proportional to the total cross-sectional area of its fibers and whose energy consumption is proportional to their volume. All of the muscle fibers in all of the motor units of a given muscle tend to move together, experiencing the same sarcomere lengths and velocities. Because of this homogeneity, most of the experimental phenomena related to contraction of whole muscle can be explained by processes occurring at the sarcomere level.

The sarcomere is the basic unit of the contractile apparatus. It is demarcated at its ends by thin Z-plates from which a matrix of thin filaments of actin project in each direction. These interdigitate with a matrix of thick filaments of myosin that are held in the center of the sarcomere by strands of the highly elastic connectin (titin) filaments that tether the thick filaments to the Z-plates. In order for muscles to contract, protrusions from myosin called myosin heads must first be cocked to a relatively high-strain configuration while bound weakly to actin and then attached more tightly to neighboring binding sites on the actin. The resulting cross-bridges act like springs that pull on the thin filaments. The metabolic energy required to cock the myosin head is provided directly by ATP molecules that are present in the sarcoplasm. In order for myosin heads to attach to neighboring binding sites, a regulatory protein called tropomyosin that normally occludes actin binding sites must undergo a conformational change. This occurs indirectly through binding of calcium to troponin, a relatively smaller molecule that is bound to tropomyosin at regular intervals along its length. When calcium binds to troponin, a local conformational change is induced that exposes nearby actin binding sites so that the cocked myosin heads can attach to form cross-bridges.

The normally relaxed state of inactive muscle is achieved by ATP-powered pumping of the calcium out of the sarcoplasm and into a network of vesicles in the sarcoplasmic reticulum called longitudinal tubules, preventing cross-bridge formation. When the muscle is activated, calcium is released from cisterns in these longitudinal tubules in response to action potentials elicited in the muscle fibers as a result of a chemical synapse with the motor axons. These action potentials propagate along the cell membrane and its invaginations deep into the muscle fiber called transverse tubules. The ATP molecules that cock the myosin heads and drive the ion pumps must be replenished eventually via oxidative catabolism of glucose and fatty acids in mitochondria, and/or glycolytic catabolism of glucose in the sarcoplasm, consuming additional energy. Some phosphate groups are hydrolyzed from the regenerated ATP and bond with creatine to form phosphocreatine (PCr) through the creatine kinase reaction. PCr functions as an additional energy reservoir for contractions and is also thought to facilitate the transport of high energy phosphates from the mitochondria to the myofilaments. The model of muscle presented here attempts to represent each of these structures and processes as explicit terms in the set of equations that comprise the model.

## Structure of the model

The two major classes of muscle models are Hill-type and Huxley models that differ mainly in their representation of the contractile element. Classical Hill-type models are phenomenological, employing arbitrary mathematical functions that relate experimental conditions (i.e. firing rate, muscle length and velocity) to the measured outcome (force) (Hill, 1938, Zajac, 1989, Thelen, 2003). Phenomenological models simply reflect the data set used to build them so they are less likely to be valid for untested conditions. Data for classical Hill models originate from maximally stimulated muscle and therefore do not account well for muscle force under physiological firing rates (Perreault et al., 2003). By contrast, Huxley-type models (Huxley, 1957, Zahalak, 1981, Ma and Zahalak, 1991) are mechanistic,; the mathematical functions that comprise them are derived directly from a hypothesized mechanism of cross-bridge dynamics. The extent to which these types of models generalize depends on the validity of the hypothesized mechanisms. Furthermore, these models tend to be complex, computationally expensive, and contain parameters that are difficult to obtain experimentally. The structure of the muscle model presented here, called Virtual Muscle, has both mechanistic and phenomenological elements to maximize generality and computational efficiency.

The model is an assembly of phenomenological models, each characterizing a major physiological process underlying muscle contraction (Figure 1). Modeling physiological processes independently as opposed to the aggregate behavior avoids over-fitting data to specific preparations and therefore improves the likelihood that the model will be valid under untested conditions. Moreover, the one-to-one correspondence of the model’s terms and coefficients to physiological processes makes it relatively straightforward to extend it to include unaccounted aspects of muscle contraction and to adjust its coefficients to model naturally occurring changes over time such as potentiation, exercise, fatigue and injury. The data driving the model originate from experiments on mammalian muscle that were designed to isolate specific processes that were then quantified. These include a wide variety of preparations with various muscle fiber architectures, but the model parameters have been normalized to dimensionless variables that depend on the constant structures of the sarcomeres, which are highly similar across all mammalian skeletal muscles. Specific models of muscle can be created by expressing their inputs in terms of these dimensionless variables (Zajac, 1989) and converting the outputs back to physical units using architectural and other muscle specific parameters (see Parameters needed to model specific muscles). The model estimates force and rate of metabolic energy consumption as a function of time in response to neural excitation, muscle length and velocity, across the full range of physiological conditions. For an example of force predictions and corresponding validation analyses, see Figures 10 and 11 in Brown and Loeb 2000; for energy predictions, see Figure 7 in Tsianos et al. 2012.

Figure 1: Model Overview. Independently modeled physiological processes and their contribution to muscle force production and energy consumption.

The model accounts for the differences among physiological processes of different fiber types, which enables simulating behavior of muscle with arbitrary fiber compositions. To model the force generating properties of individual fiber types, studies were performed on whole muscles that happen to be composed almost entirely of one fiber type. The feline soleus muscle was used for the slow-twitch model (Joyce et al., 1969, Rack and Westbury, 1969, Scott et al., 1996) and feline caudofemoralis muscle for fast-twitch (Brown et al., 1999, Brown and Loeb, 1999a,Brown and Loeb, 1999b, Brown and Loeb, 2000). The experiments were performed on anaesthetized cats at 37 degrees C with perfusion and innervation intact. Obtaining energetics data from whole muscle in vivo is challenging, so these experiments are performed typically on isolated fiber bundles. The energetics model was therefore based on in vitro studies performed on mouse soleus (slow-twitch) and extensor digitorum longus (EDL; fast-twitch) fiber bundles (Barclay, 1996, Barclay et al., 2008). These experiments were carried out at 25 degrees C because muscle fibers in these types of preparations are prone to damage at higher temperatures. Increasing the temperature has only a modest effect on energetics (Barclay et al., 2010b).

### Active force

#### Frequency-recruitment ($$U$$, $$f_{env}$$)

The drive, or synaptic current, to each motoneuron in a muscle’s motor pool is roughly the same, however, the firing rates that result in each motoneuron can be substantially different. This is because the action potentials that result in each axon depend on the change in the cell body’s membrane potential rather than the input current directly. The membrane potential induced by the input current depends on cell resistivity, which is larger for smaller motoneurons that innervate fewer muscle fibers. For this reason, as the drive increases, small motor units are recruited first, followed by larger motor units having higher recruitment thresholds (Henneman and Mendell, 1981). Increasing the drive past the threshold of a given motoneuron results in higher firing rates until the drive saturates, at which point firing rate also saturates (Monster and Chan, 1977). This monotonic relationship can be approximated by a line whose slope depends on the recruitment threshold of the motor unit (see Figure 2).

Figure 2: Motor Unit Recruitment and Modulation. The plot (top) emphasizes the relatively higher recruitment threshold of larger motor units (e.g. fast versus slow-twitch) and their firing rate modulation until the common drive (U) saturates the firing rate of all motor units. The schematic (bottom) provides a mechanistic explanation for this phenomenon, namely the size principle. The same synaptic current produces a smaller excitatory post synaptic potential (EPSP) in the larger motor neuron because its input impedance is lower. Therefore, larger input currents are necessary for overcoming the threshold for action potential generation and motor unit recruitment. Bottom figure source: Kandel ER, Schwartz JH, Jessel TM. Principles of Neural Science. 4th ed. New York (NY): McGraw-Hill; c2000. Figures 34-11; p.687.

$f_{env_{i}} = f_{min}+\left ( \frac{1-U_{MU_{i}}^{th}}{f_{max}-f_{min}} \right )*U$

Note that the firing rate is normalized to $$f_{0.5}$$ (frequency at which the motor unit produces half of its maximal isometric force) because motor units tend to fire between 1/2$$f_{0.5}$$ and 2$$f_{0.5}$$, the range over which muscle activation is most steeply modulated (see Calcium kinetics and cross-bridge activation).

#### Calcium kinetics and cross-bridge activation ($$f_{eff}$$, $$Af$$)

The steady-state relationship between motoneuron firing rate and activation (i.e. the percentage of potential cross-bridges that are attached at a given myofilament overlap) arises mainly from the release, diffusion and reuptake of calcium (see Figure 3). Calcium kinetics are modulated by other binding and buffering agents found in the muscle fiber and its sarcoplasmic reticulum, particularly parvalbumin and calsequestrin (Cannell and Allen, 1984), but these processes have not been incorporated into the model presented here. Other factors that influence the finer details of the sigmoidal relationship are outlined in Length dependency of cross-bridge dynamics. At low firing rates, the calcium released by each action potential is completely cleared before the next pulse, resulting in small, discrete twitches. As firing rate increases, the calcium released by each pulse starts to accumulate to higher concentrations depending on the rate of calcium reuptake, allowing the calcium to diffuse further and expose more cross-bridge binding sites. Eventually, the calcium concentration in all parts of the muscle fiber and throughout the entire excitation interval becomes sufficient to expose all binding sites and activation plateaus even if firing rate continues to increase. This is called a tetanic contraction.

Figure 3: Calcium Kinetics and Cross-bridge Activation. The effect of firing rate on isometric force output of muscle is shown from zero to tetanic activation. The trace is a best-fit curve for data from whole cat caudofemoralis muscle at 37 degrees C (Brown et al., 1999). Blood supply was intact and the muscle was activated via asynchronous stimulation of five bundles of motoneuron axons. Overlaid graphics highlight underlying processes at the sarcomere level that give rise to the behavior. Only a small portion of the sarcomere is shown here due to symmetry. Thick horizontal bars represent thick filaments and thin horizontal lines represent thin filaments. Vertical bars correspond to Z-disks. The small red spheres are calcium ions whose bond with actin sites is indicated by short black lines that link them and a small red dot overlaid on top of the actin (grey dots indicate inactive binding sites). If a cross-bridge is formed then the small ovals in the figure representing myosin heads are colored blue and are in a cocked configuration (large angle with respect to their neck region).

The equation below captures the effects of firing rate on cross-bridge activation (Brown et al., 1999). Y corresponds to yielding behavior exhibited by slow-twitch fibers and S corresponds to sag observed in fast-twitch fibers. Both phenomena relate to cross-bridge activation through hypothesized mechanisms discussed in the corresponding sections (see Sag and Yield). Note that equation parameter $$n_{f}$$ is length dependent (see Length dependency of calcium kinetics). The shape of the sigmoid relationship is defined by $$a_{f}$$, $$n_{f0}$$, and $$n_{f1}$$ constants, some of which are fiber type dependent. See Tsianos et al. 2012 for a complete definition of all model parameters.

$Af(f_{eff},L_{ce},V_{ce})= 1-exp\left [ -\left (\frac{YSf_{eff}}{a_{f}n_{f}} \right )^{n_{f}} \right ]$ $n_f = n_{f0} + n_{f1}\left ( \frac{1}{L_{ce}}-1 \right )$

When a motoneuron’s firing rate changes, cross-bridges attach or detach after some delay, which depends complexly on the firing rate, current level of activation and fiber length (see Length dependency of calcium kinetics). To capture this time-dependent influence of firing rate on cross-bridge activation, firing rate ($$f_{env}$$) is converted to an effective firing rate ($$f_{eff}$$) through a second order low pass filter (represented below as two first order differential equations). Time constant parameters were obtained from data and analyses performed by Brown and Loeb, 2000.

$\dot{f}_{int}\left ( t,f_{env},\bar{L}_{ce} \right )= \frac{f_{env}(t)-f_{int}(t)}{T_{f}}$ $\dot{f}_{eff}\left ( t,f_{int},\bar{L}_{ce} \right )= \frac{f_{int}(t)-f_{eff}(t)}{T_{f}}$ $T_{f}= \left\{\begin{matrix} T_{f1}\bar{L}_{ce}^{2}+T_{f2}f_{env}(t),\dot{f}_{eff}\geq 0 \\ \left (T_{f3}+T_{f4}Af \right )/\bar{L}_{ce},\dot{f}_{eff}< 0 \end{matrix}\right.$

#### Sag ($$S$$)

When fast-twitch muscle is excited isometrically with a constant drive, its force output increases initially to some maximal level and then decreases gradually (Burke et al., 1973, Brown and Loeb, 2000; Figure 4). This is thought to occur due to an increase in the rate of calcium reuptake, which effectively reduces the concentration of calcium in the sarcoplasm and hence the number of crossbridges that can form (see Relevant Physiology).

Figure 4: Sag. Experimental data depicting the sag phenomenon in response to constant stimulation for several frequencies. The traces depict actual data recorded from whole cat caudofemoralis muscle at 37 degrees C (Brown et al., 1999). Blood supply was intact and the muscle was activated via asynchronous stimulation of five bundles of motoneuron axons. Force for each stimulation frequency is normalized to the force generated at 133ms, which was used to define the $$A_{f}$$ relationship (see Calcium kinetics and cross-bridge activation). The graphics on the bottom portion illustrate the effect of increasing the rate of calcium reuptake on calcium concentration and cross-bridge formation.

This phenomenon can be modeled by a time-varying scaling factor applied to $$f_{eff}$$, which is directly related to the calcium present in the sarcoplasm (see Calcium kinetics and cross-bridge activation; Brown and Loeb, 2000).

$\dot{S}(t,f_{eff}) = \frac{a_s-S(t)}{T_s}$

#### Force-length ($$FL$$)

The amount of active force a muscle produces during maximal isometric contractions (i.e. muscle length held constant during the contraction) depends on the length at which it is fixed (Figure 5). It produces large forces at intermediate lengths and relatively smaller forces at either shorter or longer lengths (Gordon et al., 1966, Scott et al., 1996).

At some intermediate length, typically referred to as optimal length, overlap between actin and myosin filaments is maximal. Therefore, the number of cross-bridges that could form and the contractile force are also maximal. Myofilament overlap is less at longer lengths so the number of potential cross-bridges and force are also smaller. Muscle force decreases at incrementally longer lengths until the point at which no more cross-bridges can form and muscle can no longer produce active force. At relatively shorter lengths, actin filaments start to slide passed each other (i.e. double overlap of actin), which is thought to sterically hinder the formation of cross-bridges, therefore, reducing the amount of force that can be produced. At lengths less than 0.75 $$L_{0}$$, the thick filament is thought to collide with the Z-disks, thereby generating a passive restoring force opposing the contractile force (Brown et al., 1996b, Brown et al., 1999). This passive component of the force-length relationship is modeled separately and described in Thick filament compression. At lengths beyond $$L_{0}$$, muscle connective tissue is pulled taut and produces a restoring force in the direction of the contractile force. This passive component is also modeled separately and described in Parallel elastic element.

Figure 5: Force-Length relationship. Muscle force dependence on fascicle length is shown for tetanically stimulated muscle. The trace is a best-fit curve for data from whole cat soleus muscle at 37 degrees C (Brown et al., 1996b, Scott et al., 1996) with the effects of thick filament compression removed (see text for details). Nerve and blood supply was intact and the muscle was activated via cuff electrodes placed around the tibial nerve. Overlaid graphics depict myofilament overlap and effects on cross-bridge formation for different lengths. Note that the FL relationship does not include the steeper decline at the shortest sarcomere lengths that arises from thick-filament compression, as described below.

The force corresponding to optimal length is by definition maximal so normalizing all force profiles to maximal isometric force under tetanic conditions provides a useful reference point. The equation capturing the normalized relationship is shown below.

$FL(L_{ce}) = exp\left ( -\left |\frac{{L_{ce}}^\beta - 1}{\omega } \right |^{\rho }\right )$

ω, β, ρ are constants that define the exact shape of the inverted-U and are all fiber-type dependent.

#### Force-velocity ($$FV$$)

Maximally stimulated muscle at a particular length produces different levels of force depending on the rate at which it is shortening or lengthening; it produces less force while shortening (Hill, 1938) and more force while lengthening, relative to the isometric condition (Gasser and Hill, 1924, Scott et al., 1996; Figure 6).

Shortening muscle is accompanied by motion of the attached cross-bridges toward relatively less strained configurations (Ford et al., 1985). Larger rates of shortening increase the probability that a given cross-bridge will be in a low-strain configuration; therefore, the force generated by the entire population is lower overall. By contrast, cross-bridges in lengthening muscle become relatively more stretched (Piazzesi et al., 1992) and generate more tension until they are forcibly ripped away from the actin binding sites.

Figure 6: Force-Velocity relationship. Dependence on muscle force output on the velocity of fascicle stretch (V>0) and shortening (V<0) for tetanic stimulation at optimal muscle length. The trace is a best-fit curve for data from whole cat soleus muscle at 37 degrees C (Brown et al., 1996b, Scott et al., 1996). Nerve and blood supply was intact and the muscle was activated via cuff electrodes placed around the tibial nerve. The graphics illustrate the hypothesized mechanisms for force enhancement relative to isometric in the lengthening case and depression in the shortening case. Note that in reality myosin heads do not move synchronously as shown in the figure because they are subject to thermodynamic noise and other stochastic events; the graphic illustrates their configuration on average.

The equation below is included in the model to account for this phenomenon.

$FV(V_{ce},L_{ce})=\left \{ \begin{matrix} (V_{max}-V_{ce})/[V_{max}+(c_{v0}+c_{v1}L_{ce})V_{ce}],\; \; \; \; \; \; \; V_{ce}\leqslant 0\\ [b_v-(a_{v0}+a_{v1}L_{ce}+a_{v2}{L_{ce}}^{2})V_{ce}]/(b_v+V_{ce}),V_{ce}>0 \end{matrix} \right.$

$$c_{v0}$$, $$c_{v1}$$, and $$V_{max}$$ are constants that define the shortening end of the relationship and depend on fiber type. $$V_{max}$$ corresponds to the maximum velocity of shortening that the muscle fiber can undergo. The lengthening portion of the relationship is defined by fiber type dependent constants $$b_{v}$$, $$a_{v0}$$, $$a_{v1}$$, and $$a_{v2}$$. Note that the force-velocity relationship also depends on muscle length. See Potentiation for a detailed description of this phenomenon as well as a mechanistic explanation.

The difference in tetanic force relative to isometric for shortening and lengthening conditions may also be attributed to some extent to a change in the number of crossbridges attached. In fact, it has recently been concluded that the force-velocity relationship can be explained almost entirely by a change in cross-bridge number rather than cross-bridge strain (Brunello et al., 2007, Piazzesi et al., 2007, Barclay et al., 2010a, Fusi et al., 2010). This interpretation, however, is based on estimates of fiber stiffness that may be inaccurate. Muscle fiber stiffness is measured by applying length perturbations to activated muscle (Ford et al., 1977) that are small (less than the critical shortening and lengthening displacements that would cause a significant number of crossbridges to detach) and rapid (faster than the time it would take for crossbridges to actively change their configuration). Stiffness is then simply calculated as the magnitude of the force change divided by the magnitude of the imposed length change. To estimate the contribution of crossbridge stiffness, which is proportional to the number of crossbridges attached, it is often assumed that myofilament stiffness can be modeled as a linear spring in series with another spring representing the lumped influence of all crossbridges on stiffness. However, it is likely that myofilament stiffness is nonlinear (Irving et al., 2011), which has been shown to have a large influence on estimates of number of attached crossbridges (Mansson, 2010). Crossbridge stiffness has also been shown to be highly nonlinear from measurements performed at both the molecular (Kaya and Higuchi, 2010) and fiber level (Nocella et al., 2013). Nocella et al. performed high resolution measures of fiber stiffness during the types of muscle length perturbations described above and observed that stiffness changed with force synchronously during muscle fiber length perturbations. It is plausible that stiffness changes during shortening and lengthening contractions reflect a change in cross-bridge stiffness resulting from change in strain rather than a change in the number of cross-bridges attached.

It is also worth noting that these interpretations of the Force-Velocity relationship are based on experiments on frog muscle held at a temperature near 0 degrees C, which is substantially lower than physiological temperature. The Force-Velocity relation is affected significantly by temperature (Rall and Woledge, 1990). A given increase in shortening velocity leads to a smaller decrease in force for higher temperatures. This could reflect a smaller effect of contractile velocity on the number of crossbridges attached at this temperature, which can be explained by an increase in the rate of crossbridge attachment. It may also reflect cross-bridges being in a relatively stiffer configuration, given the larger strain experienced by the cross-bridges at higher temperatures (Ford et al., 1981, Colombini et al., 2008).

In either case, the underlying mechanism will not affect the validity of the force model because the force-velocity relation is obtained directly from empirical data rather than derived from a hypothesized mechanism. The number of crossbridges formed does affect the energetics of contraction, however. For submaximal levels of activation, slow-twitch muscle could experience yielding where the number of cross-bridges decreases substantially for both the shortening and lengthening condition. This phenomenon and model have been described in Yield.

#### Interactions among length, velocity and activation

##### Length dependency of cross-bridge dynamics

Cross-bridge activation increases with length. Theoretically, this can result from either an increase in the number of cross-bridges attached or an increase in cross-bridge strain, but cross-bridge strain has been shown to be similar across sarcomere lengths (Gordon et al., 1966). The number of cross-bridges increases with higher rate of cross-bridge attachment, lower rate of cross-bridge detachment, or higher concentration of calcium in the sarcoplasm. Calcium concentration, however, has been shown to be independent of length (Balnave and Allen, 1996).

Changing length of the sarcomeres has effects on the rate of attachment and detachment of cross-bridges to activated and exposed binding sites on the actin. Because a muscle fiber does not gain or lose volume as it changes length, the fiber must have a larger diameter when it is at a shorter length. The hexagonally packed lattice of myofilaments in each sarcomere will then be more widely spaced, changing the distance between the myosin heads and the thin filaments where they must bind to form cross-bridges. The heads are located on the ends of hinged arms containing myosin light chains, which are canted away from the longitudinal axis of the thick filament. When the sarcomere is at a long length, the myosin heads barely fit between the thick and thin filaments, so they are close to and can bind rapidly to form cross-bridges. The small distance between the myosin heads and actin binding sites also increases the strength of the actomyosin bond, which would reduce the probability of cross-bridge detachment. At short lengths, the myosin heads are less favorably disposed and their binding may also be affected adversely by the double-overlap of the thin filaments. These effects are particularly large at lower levels of activation, where activated binding sites are more scarce on the thin filaments. The net result is shown in Figure 7, which is based on data and models first described in Brown et al., 1999.

Figure 7: Length Dependence of Cross-bridge Activation. Cross-bridge activation is measured as the percentage of maximal force generated for a given amount of myofilament overlap (i.e. effects of myofilament overlap are removed from muscle force). The trace is a best-fit curve for data from whole cat caudofemoralis muscle at 37 degrees C (Brown et al., 1999). Blood supply was intact and the muscle was activated via asynchronous stimulation of five bundles of motoneuron axons. Fascicle length was fixed at 0.8, 1, and 1,2L0. In general, a larger portion of the available contractile machinery is activated for the same firing rate when the muscle is stretched.

The precise slope of this relationship, i.e. the sensitivity of force to firing rate, is affected by additional factors. A low firing rate activates only a small portion of the binding sites on the actin filaments. An incremental increase in the firing rate results in a relatively larger increase in the number of cross-bridges formed. This is probably due to a relatively larger amount of exposed binding sites per active troponin molecule; this mechanism is known as cooperativity (Gordon et al., 2000, Loeb et al., 2002). Activation of troponin induces a local conformational change in tropomyosin, exposing neighboring binding sites on the actin. When more troponin molecules are activated along tropomyosin, a more global conformational change could be occurring, thus, freeing up more actin binding sites. Crossbridge formation itself may also facilitate exposure of adjacent binding sites by inducing relative motion between actin and tropomyosin.

##### Length dependency of calcium kinetics

Changing length of the sarcomeres has effects on the calcium kinetics that govern activation. The cisterns from which the calcium is released appear to be tethered to the Z-plates at a location that is near the middle of the actin-myosin overlap when the muscle is at optimal length (Brown et al., 1998). At longer lengths, calcium has to diffuse over longer distances to get to the actin binding sites, which takes more time. Thus the rise times of the activation are longer (see Brown et al., 1999, Brown and Loeb, 2000). As mentioned in the previous section (Length dependency of cross-bridge dynamics), it is likely that cross-bridge attachment rate increases and detachment rate decreases with length and would contribute to the decrease in muscle relaxation time observed at longer lengths (Brown and Loeb, 1999a, Rassier and Macintosh, 2002).

##### Potentiation

It was originally observed that the force generated by a single twitch of a fast-twitch muscle was much larger after a brief tetanic contraction than before (Hughes, 1958). In fact, the potentiated state tends to prevail for minutes after a few, brief trains of stimulation at physiological rates, suggesting that the normal operating state of fast-twitch muscle is the potentiated state and that the conditions observed after a long period of quiescence reflect a dispotentiated state (Brown and Loeb, 1998). Furthermore, the dispotentiated state exhibits rather odd behavior such as a pronounced rightward shift in the shape of the force-length curve at subtetanic frequencies. For these reasons, the model of the fast-twitch muscle fibers presented here captures their behavior in the fully potentiated state.

Potentiation and dispotentiation phenomena appear to be related mainly to the effects of myofilament lattice spacing described above, but see also (Smith et al., 2013). During dispotentiation, the light chain becomes dephosphorylated and the cant angle becomes smaller, leaving the myosin heads further from the thin filaments (see Stull et al., 2011 for review). The effect is greater at short lengths when the lattice spacing is large. This adversely affects the kinetics of cross-bridge formation in general, but more so at low firing rates when exposed binding sites are more sparsely distributed. After a few cycles of calcium release in an active muscle, the short chains become fully phosphorylated and the cant angle increases so as to position the myosin heads more favorably at short and optimal muscle lengths. At longer lengths and narrower lattice spacings, the myosin heads probably remain favorably positioned because they cannot get by the closely spaced thin filaments. The force-length relationship thus returns to the simple inverted-U with maximum at optimal myofilament overlap regardless of firing rate. Slow twitch muscle appears to lack such a dephosphorylation process and so behaves always as if fully potentiated. The function of the dispotentiation remains obscure but may be related to minimizing “stiction” (Matthews, 1981), which arises when stray cross-bridges in an inactive muscle form spontaneously and resist passive stretching by antagonist muscles.

The effects of sarcomere length on myofilament lattice spacing give rise to length-dependencies of the force-velocity relationship that are different for slow- and fast-twitch muscle (Figure 8). For slow twitch muscle, the force-velocity relationship depends on length only for the lengthening condition (Scott et al., 1996), while fast twitch muscle exhibits a dependency for the shortening condition (Brown and Loeb, 1999a). In lengthening slow-twitch muscle, smaller spacing between myosin heads and binding sites would result in a higher probability of crossbridge formation, thus, leading to higher forces produced for the same velocity of stretch. This effect is also observed in dispotentiated fast-twitch muscle (see Potentiation), but not in the potentiated state, because the reduction of interfilament spacing resulting from potentiation is thought to mask the effects of length changes (Brown and Loeb, 1999a). In shortening fast-twitch muscle, it is thought that smaller interfilament spacing at longer lengths increases the probability that a given crossbridge remains attached, thus, making it more likely that it occupies relatively low-strain, low force configurations. Crossbridges in slow-twitch muscle stay attached for longer periods of time, which may mask the effects of length on the rate of detachment hypothesized for fast twitch muscle.

Figure 8: Length Dependence of Force-Velocity. Force-Velocity curves depicting a substantial fascicle length dependency of slow-twitch muscle in the lengthening state (left) and fast-twitch muscle in the shortening state (right). Traces for L = 0.8, 1, and 1.2 L0 are shown. The traces in the slow-twitch plot are best-fit curves for data from whole cat soleus muscle at 37o C (Brown et al., 1996b, Scott et al., 1996). Nerve and blood supply was intact and the muscle was activated via cuff electrodes placed around the tibial nerve. The traces in the fast-twitch plot are best-fit curves for data from whole cat caudofemoralis muscle at 37 degrees C (Brown et al., 1999). Blood supply was intact and the muscle was activated via asynchronous stimulation of five bundles of motoneuron axons.

##### Yield ($$Y$$)

When slow twitch muscle is submaximally stimulated at constant length and then stretched or shortened, its force output declines (Joyce et al., 1969). This effect intensifies for lower stimulation frequencies, as shown in Figure 9. The mechanism for this phenomenon has been hypothesized to be a reduction in the number of crossbridges attached. As muscle length changes past the stroke length of the attached crossbridges, myosin heads detach from their binding sites and must reattach to continue to produce force. The rate at which slow-twitch myosin heads can attach appears be relatively slow compared to the rate at which binding sites slide past them, therefore, making it difficult for crossbridges to reform and produce force. At lower frequencies of stimulation, relatively few binding sites are available on the helical actin, making it less likely that they will be in an appropriate orientation and distance from myosin heads for binding. By applying a time-varying scaling factor to the effective firing rate parameter (see Calcium kinetics and cross-bridge activation), the model is able to capture this effect well (Figure 9).

$\dot{Y}(t)=\frac{1-c_Y\left [ 1-exp\left ( -\frac{\left |V_{ce} \right |}{V_Y} \right ) \right ]-Y(t)}{T_Y}$

$$c_{Y}$$ and $$V_{Y}$$ are constants that define the intensity of the yielding effect and $$T_{Y}$$ defines its time course.

Figure 9: Yielding. Plot of Force-Velocity relationships of slow-twitch muscle for various frequencies of stimulation. As opposed to tetanic stimulation, the force-velocity relationships for submaximal firing rates are marked by a sharper decline in force in both shortening and lengthening states relative to isometric. Experimental data and model predictions, with and without accounting for yielding, are shown. Experimental data originated from whole soleus muscle at normal body temperature (Joyce et al., 1969). Figure source: Fig 10c from Brown et al. 1999.

### Passive elastic elements

Muscle is comprised of many elastic components that generate passive forces in response to tensile or compressive strain. The aggregate force of these components may be substantial relative to the active force depending on the muscle’s condition of use.

#### Parallel elastic element ($$F_{pe1}$$)

Inactive muscle is under tension for about half of its anatomical range of lengths. Its passive tensile force rises exponentially over a relatively low range of lengths followed by a steeper linear increase that persists all the way until maximum anatomical length (Figure 10). The precise nature of the relationship for a particular muscle depends mainly on its maximum isometric force and maximum anatomical length. It is negligible over most of the range and peaks at less than 10% of $$F_{0}$$ at maximum anatomical length.

Figure 10: Parallel Elastic Force. Passive muscle force is plotted for a representative muscle having a maximum anatomical length of 1.1L0. This relationship was derived from passive force measurements performed on five strap-like muscles of the cat (Brown et al., 1996a). In general, muscles with larger anatomical ranges of motion produce less passive force at the same length (relative to optimal), that is, their passive force-length curve shifts to the right. As depicted by the superimposed schematics, some of this force results from stretching an elastic element in the sarcomere called connectin (also titin) that links the thick filament to the Z-disks, but see text for role of endomysial connective tissue. Myosin heads and actin binding sites are omitted for clarity.

In frogs, this passive tension appears to arise primarily from stretching of the spring-like connectin myofilaments that tether the myosin molecules to the Z-disks (Magid and Law, 1985). The contribution of intermediate filaments like desmin that interconnect myofibrils within a muscle fiber has been shown to be negligible over the physiological range of sarcomere lengths (Wang et al., 1993). Compared to frogs, mammalian muscles have variable but considerably more endomysial collagen surrounding each muscle fiber. This may account for substantial variation in their passive force curves, which can be substantially shifted to the left or right with respect to the optimal sarcomere length defined by myofilament overlap (Brown et al., 1996a). This requires introduction of an additional morphometric parameter into the muscle model, $$L_{max}$$, which is based on the maximal anatomical length experienced by the muscle when attached to the skeleton. Below is the equation capturing the effect of this property on the entire passive force-length relationship.

$F_{pe1}(L_{ce},V_{ce})=c_1k_1ln\left \{ exp\left [ \frac{L_{ce}/L_{ce}^{max}-L_{r1}}{k_1} \right ]+1 \right \}+\eta V_{ce}$

$$c_{1}$$, $$k_{1}$$, and $$L_{r1}$$ are constants defining the precise shape of the curve. Linear damping was also incorporated into the equation with constant coefficient, η, to account for the small amount of damping observed experimentally, which also contributes to the mathematical stability of the model.

#### Thick filament compression ($$F_{pe2}$$)

At very small lengths, muscles generate passive force that opposes the force of contraction. This force is negligible for long lengths but starts to increase exponentially from 0.7$$L_{0}$$ until the minimum physiological length where net contractile force goes to zero (see Figure 11; Brown et al., 1996b).

The steady-state force produced by tetanically stimulated muscle decreases at a fairly constant rate for the range of lengths from $$L_{0}$$ and 0.7$$L_{0}$$, reflecting steric interference between the overlapped thin filaments. At this relatively small length, however, the slope steepens abruptly and remains steep until the length at which net muscle force is zero. Assuming that the effects of double-overlap on cross-bridge formation are the same over the entire range of lengths in which it occurs, the number of cross-bridges, hence force, should decrease proportionally with decreasing length. The abrupt change in force appears to reflect a collision between the myosin filament and the Z-disk, which would generate a reaction force opposing the force of contraction. The extent of compression of the myosin filament onto the Z-disk would be less for submaximal contractions because the number of cross-bridges generating force would be relatively smaller. Thus it is important to model this effect separately from the force-length relationship, which affects force generation equally at all activation levels. The equation below computes the force at the muscle level as a function of length that opposes the force of contraction.

$F_{pe2}(L_{ce})=c_2\left \{ exp[k_2(L_{ce}-L_{r2})]-1 \right \}$

Constants $$c_{2}$$, $$k_{2}$$, and $$L_{r2}$$ define the exact shape of the relationship depicted in Figure 11.

Figure 11: Thick Filament Compression. Plot shows the parallel passive force of muscle opposing the contractile force as a function of fascicle length. The schematics illustrate the hypothesis that this force arises due to compression of the thick filament as the Z-disks are pulled even closer together at short muscle lengths. The precise form of this relationship was derived from force-length measurements of whole cat soleus muscle at 37o C (Brown et al., 1996b, Scott et al., 1996). Nerve and blood supply was intact and the muscle was activated via cuff electrodes placed around the tibial nerve.

#### Series elastic element (tendon + aponeurosis; $$F_{se}$$)

Muscles exert forces on bone segments through intermediate connective tissue, comprised mostly of aponeurosis and tendon. This serially connected tissue affects the muscle’s motion, which modulates force output substantially. Clearly, the length of a muscle for a particular joint configuration depends on the length of the series elastic element because together they must span the entire path from the origin to the insertion site of the musculotendon (see Musculoskeletal Mechanics and Modeling). The portion of the path that is occupied by muscle depends not only on the resting length of the series elastic element, but also on its elasticity, which determines how much it will stretch when it is pulled by the attached muscle. If a muscle is activated and then deactivated under isometric conditions as defined at the origin and insertion, the sarcomeres may experience large changes in length and especially velocity that will affect their force generating ability. Under dynamic conditions of activation and kinematics, the muscle fibers may actually move out of phase with the musculotendon (Zajac et al., 1981). This has important implications for the ability of muscle spindles to encode musculotendon length, hence joint position (Hoffer et al., 1989; see Proprioceptors and Models of Transduction).

The magnitude of the effects described above depends on the ratio of the lengths of the muscle fascicles to the connective tissue in series with them. Highly pinnate muscles have short muscle fibers oriented obliquely between fascial planes called aponeuroses that merge with the external tendon. Proximally, the aponeurosis has a thin sheet-like structure with a large surface area to accommodate insertion of a large number of fibers while more distally its shape gradually becomes thicker and smaller in surface area. It has been shown experimentally that for a wide range of muscle forces, aponeurosis and tendon strain are the same (Scott and Loeb, 1995). The fact that the gradual thickening of the aponeurosis is paralleled by an increase in the number of fibers that exert force on it suggests that the stresses are fairly constant along its entire length. Because the stress and strain experienced by the aponeurosis and tendon are roughly the same, then presumably their material properties are similar and they can be lumped into a single element. The absolute thickness and strength of this connective tissue element is presumably controlled by trophic factors that depend on the maximal tensile forces that it experiences. The model assumes that tendon force scales with muscle cross-sectional area, or equivalently it scales with $$F_{0}$$, and strain is measured with respect to tendon length at $$F_{0}$$ ($$L_{0}^T$$), permitting a generic relationship to be derived (see Figure 12). This has been confirmed using several force-length relationships of the tendon/aponeurosis obtained experimentally (Brown et al., 1996b). $$L_{0}^T$$ is a good normalization factor because it can be determined more reliably as opposed to the more commonly used tendon slack length. The equation that defines this relationship is shown below.

$F_{se}(L_{se})=c^Tk^Tln\left \{ exp\left [ \frac{L_{se}-L_{r}^{T}}{k^T} \right ]+1 \right \}$

Constants $$c_{T}$$, $$K_{T}$$, and $$L_{r}^{T}$$ define the precise form of the relationship shown in Figure 12.

In reality, the entire aponeurosis is not in series with any given muscle fiber, because each fiber exerts force at a distinct location along its length. The force exerted by a given fiber on its attachment site on the aponeurosis acts in parallel with the forces transmitted via other parts of the aponeurosis. It has been shown that this arrangement can yield substantially different musculotendon behavior theoretically (Epstein et al. 2006), but such behavior has not been shown to occur in vivo (Tilp et al. 2012).

Series elastic compliance of the myofilaments of the muscle fiber has a small effect on sarcomere length during isometric contractions (Kawakami and Lieber, 2000); therefore, it does not affect force generation through changes in myofilament overlap. Theoretically, relative motion of the myofilaments may occur at the nm scale, which would affect the energetics through an increase in cross-bridge cycling rate (see Energy related to cross-bridges). The effects on force output are likely to be small because while fluctuations in sarcomere length may cause some cross-bridges to detach, it may bring detached myosin heads to a favorable position for binding to actin sites. Moreover, while the shortening phase of the oscillation would reduce the strain on the cross-bridges and therefore reduce the force they generate, the lengthening phase would increase strain, hence their force output. In physiological contractions, cross-bridges are formed asynchronously throughout the muscle so force and strain changes are likely to be balanced, leading to a small net change in force attributed to series elastic compliance.

Figure 12: Series Elastic Tendon/Aponeurosis Force. Force generated by elastic tissue that is in series with muscle as a function of its length. The relationship was derived from data from whole soleus muscle of the cat (Scott and Loeb, 1995). The schematic shows that for a given musculotendon length and low muscle activation (bottom), the tendon will be relatively slack and therefore produce little restoring force. At higher activation levels (top), the muscle pulls more on the tendon and causes a larger restoring force. The stiffness is also greater in this case, denoted by the increased slope of the relationship, because of the various conformational changes in the collagen that obtain at different strains (Diamant et al., 1972).

### Energetics

Muscle contraction is an active process whose energy is supplied directly by the high-energy phosphate bonds of ATP molecules ($$E_{initial}$$). The expended ATP molecules must then be replaced via metabolism in the muscle to maintain fuel supply for future contractions; this process requires additional energy ($$E_{recovery}$$). Researchers measure the rate at which muscles consume energy experimentally by summing the heat output rate and mechanical power associated with brief contractions under various conditions. Characterizing the rate of change of energy consumption under each condition independently enables the prediction of total energy loss for situations in which all conditions are changing concurrently.

#### Energy driving the contraction ($$E_{initial}$$)

To measure the energy that activates the contractile machinery directly, researchers typically investigate relatively brief contractions that have a shorter duration than the time it takes for metabolic pathways to begin consuming significant amounts of energy. Most of this energy fuels two major physiological processes that underlie contraction: calcium reuptake into the sarcoplasmic reticulum ($$E_{a}$$) and cross-bridge cycling ($$E_{xb}$$).

#### Energy related to ion pumps ($$E_{a}$$)

When a muscle fiber is excited at a particular firing rate, voltage gated channels in the sarcoplasmic reticulum release calcium ions down their concentration gradient into the sarcoplasm, where they activate the contractile machinery. Eventually the muscle will relax, which requires all of the calcium to be resequestered to deactivate the cross-bridges and restore the original concentration in the sarcoplasmic reticulum. Transporting calcium ions against their concentration gradient is an active process that is mediated by ATP activated pumps that are distributed all along the sarcoplasmic reticulum’s longitudinal tubules. This energy ($$E_{a}$$) can be isolated in thermodynamic experiments by measuring $$E_{initial}$$ when cross-bridges are not allowed to form. This can be accomplished pharmacologically (using N-benzyl-p-toluenesulphonamide) or by stretching the muscle to the point where myofilaments do not overlap. The energy consumed by this process is mainly a function of firing rate, independent of sarcomere length and velocity. It has been shown experimentally that Ea is roughly one third of the energy related to cross-bridge cycling when muscle is fixed at its optimal length (see Figure 13). Although incomplete, evidence suggests that this is true over a large range of firing rates, which means that the $$E_{a}$$ vs. firing rate relationship ought to have a similar shape to cross-bridge activation vs. firing rate in the isometric case. Thus, the $$E_{a}$$-firing rate relationship is modeled using the cross-bridge activation function that is scaled to convert it to units of energy (W). The appropriate scaling factor is defined by constants $$e_{3}$$ and $$e_{4}$$ in the equation below. Given its sigmoidal relationship, this component of energy consumption is most sensitive over an intermediate range of firing rates and least sensitive at either extreme.

$E_a(f_{eff}) = \frac{e_3}{3e_4}Af(f_{eff},L_{ce}=L_{ce0},V_{ce}=0)$

#### Energy related to cross-bridges ($$E_{xb}$$)

The pulling force generated at each end of a muscle represents the sum of many smaller pulling forces generated by a vast number of constituent contractile elements (sarcomeres and their myofilaments) arranged both in series and in parallel. During a contraction, pulling forces are generated on both ends of the sarcomeres synchronously that are transmitted to the ends of the muscle through connective tissue. If the external load is lower than the contractile force, then the Z-disks will come closer together and the muscle will shorten and perform mechanical work on the load, and the sarcomeres will perform work on the connective tissue. The myofilaments will slide past each other and the cocked myosin heads will dissipate their stored mechanical energy as they perform this work. When a cross-bridge reaches the angle at which it exerts little or negative force on the thin filament, the myosin head simultaneously binds a new molecule of ATP, desphosphorylates it to extract its chemical energy, which enables it to detach from the thin filament and recock. If the external load is greater than the contractile force, the cross-bridges are stretched until they are forcibly detached, consuming no additional ATP and remaining in the cocked state, so they can immediately reattach to the next available binding site. The energy consumption related to cross-bridge cycling is thus highly dependent on velocity of sarcomere motion (see Figure 13). It increases rapidly with increasing velocity of shortening (i.e. negative velocity) and it falls rapidly to zero for positive velocity of stretch. Under isometric conditions (zero velocity), cross-bridges continue to cycle at a non-zero rate because of the squirming motion of the myofilaments caused by thermodynamic noise and the compliant nature of the cross-bridges, actin, and other structural proteins in series.

Figure 13: Partitioning of Energy Driving Muscle Contraction. Total rate of energy consumed by muscle (Einitial) as a function of its velocity as well as its partitioning into energy related to excitation $$E_{a}$$ and energy related to cross-bridge cycling $$E_{xb}$$ components is shown for slow and fast-twitch fibers. These relationships were derived mostly from small bundles of mouse muscle fibers (Soleus for slow-twitch and EDL for fast-twitch) at 25 degrees C (Barclay, 1996, Barclay et al., 2008). For the most part, muscles were stimulated maximally at their optimal length. See Tsianos et al. 2012 for a complete derivation of these relationships.

Energy related to cross-bridge cycling also depends on length and firing rate, which together determine the portion of the contractile machinery that is active, i.e. the number of cross-bridges that are cycling. The cross-bridge activation (Af) and force-length (FL) relationships together estimate the percentage of potential cross-bridges that form so they can be used to scale the $$E_{xb}$$-velocity relationship obtained at tetanic stimulation and optimal length, i.e. the condition at which all potential cross-bridges are formed (Af*FL = 1).

The $$E_{xb}$$ relationship at tetanic stimulation and optimal length ($$E_{xb}$$ curve in Figure 13) is defined by the piece-wise function in the equation below. The relationship for the shortening condition is approximated well by a rational function while the lengthening condition is approximated by a line. Muscles cannot convert mechanical energy to chemical energy (i.e. metabolic energy consumption cannot be negative); therefore, $$E_{xb}$$ is equal to zero when lengthening velocity exceeds a critical value defined as $$V_{ce0}$$. Constants $$e_{1}$$, $$e_{2}$$, $$e_{3}$$, and $$e_{4}$$ define the precise form of the relationship. Note that $$e_{3}$$ and $$e_{4}$$ are the same constants used in the equation of $$E_{a}$$.

$E_{xb}(f_{eff},L_{ce},V_{ce})=Af(f_{eff},L_{ce},V_{ce})FL(L_{ce})*\left\{\begin{matrix} \frac{e_1{V_{ce}}^{2} + e_2V_{ce}+e_3}{e_4-V_{ce}}-\frac{e_3}{3e_4},\; \; V_{ce}\leq 0\\ \frac{2e_3}{3e_4} +\left ( \frac{e_2e_4+e_3}{{e_4}^2} \right )V_{ce}, \; \; 0<V_{ce}\leqslant V_{ce0} \\ 0, \; \; V_{ce}>V_{ce0} \end{matrix}\right.$

$where\; \; V_{ce0}=-\left ( \frac{2e_3}{3e_4} \right )\left ( \frac{{e_{4}}^{2}}{e_2e_4+e_3} \right )$

#### Energy related to metabolism ($$E_{recovery}$$)

The ATP used for activating the contractile machinery must be recovered to maintain the fuel supply for subsequent contractions. This is accomplished via metabolism within the muscle, which lags the contraction by roughly one second and lasts for about two minutes following relaxation. Energy related to metabolism is quantified by eliciting a muscle contraction over a sufficiently brief period (less than one second) and measuring all of the energy consumed after the muscle relaxes. $$E_{recovery}$$ is typically expressed as a function of $$E_{initial}$$ in the literature because the level of metabolic activity naturally depends on the fuel supply that needs to be restored.

The $$E_{recovery}$$/$$E_{initial}$$ ratio has been investigated mainly for tetanic contractions, but it is reasonable to assume that it is the same for submaximal contractions because theoretically both $$E_{recovery}$$ and $$E_{initial}$$ should be directly proportional to the number of ATP molecules expended, assuming the metabolic pathways and cellular environment are identical. The ratio has been shown to differ substantially between slow and fast-twitch fibers (1.5 versus 1, respectively) undergoing oxidative metabolism. The ratio for glycolytic metabolism has not been investigated rigorously, but theoretical predictions suggest that it is similar. These predictions are based on the energy wasted as heat per ATP molecule during catabolism of glucose. Glycolytic metabolism, however, entails the additional energetic cost of transporting the lactate byproduct to the liver, converting it to pyruvate, and redistributing for subsequent catabolism by oxidative pathways. This cost is likely to be substantial but has not been quantified yet.

Because metabolism occurs more slowly than the contraction itself, its dynamics must be modeled in order to account for the amount of energy consumed at any given time as opposed to the total energetic cost of a contraction. Leijendekker and Elzigna (1990) investigated the rate of energy consumed over time by slow and fast-twitch muscle following brief contractions. Assuming that the dynamics of metabolism are linear, the rate of recovery energy consumed over time can be estimated for an arbitrary contraction by convolving the energy driving the contraction ($$E_{initial}$$) with the system’s response to an $$E_{initial}$$ impulse, which can be derived easily from the experiment (see Tsianos et al., 2012). Using this approximation, the model does a remarkable job at predicting the dynamics of recovery energy for human muscle contracting in vivo.

## Balancing accuracy and computational efficiency

The feasibility of studying large-scale neuromusculoskeletal systems using models depends largely on their accuracy and computational efficiency. The model presented here is highly accurate, but because it models the dynamics of each motor unit independently, it requires many computations to solve for all of their states at any given time. To reduce the model’s computational load, an algorithm was devised that mathematically lumped all motor units of each fiber type into one. This improved the model’s efficiency approximately 10-fold without having significant effects on its accuracy (Tsianos et al., 2012).

The lumped model computes an effective firing rate that would reproduce the aggregate force output of a realistic number of motor units with a range of firing rates associated with any given excitation level of the motor pool. The weighting function used to compute the firing rate of the lumped unit is dependent on the drive to the muscle because different levels of drive cause motor units to fire at different levels and may also recruit a different set of motor units (see Frequency-recruitment). The dynamics of force generation also depend on the number of motor units firing and their relative activation as well as each of their activation histories. The lumped unit’s dynamics must therefore be modulated by the level of drive to capture this effect. It turns out that a first-order low pass filter applied to the drive with a time constant that is a function of the change in drive is sufficient. It captures most of the complex dynamics and corresponding force output of the validated muscle model for a wide range of excitation conditions (see Figure 10 in Tsianos et al., 2012). This modification to the model improved the computational efficiency substantially because modeling calcium kinetics originally required solving two states per motor unit. A realistic model of muscle having 50 motor units would have 100 states (2*50) while the lumped model has 3 (2+1 related to the additional first order filter applied to the drive signal).

## Parameters needed to model specific muscles

Musculotendons are comprised of the same basic elements, such as sarcomeres for contractile elements and collagen for passive elements, whose structural organization and properties lead to substantially different behavior on the whole musculotendon level. Knowledge of musculotendon morphometry and fiber-type specific properties of sarcomeres is crucial for tailoring the generic relationships that comprise the model to the specific muscle of interest.

### Optimal fascicle length ($$L_{0}$$)

As mentioned in Relevant physiology, sarcomeres have roughly the same dimensions and experience the same range of length changes during a contraction. Longer muscle, therefore, has more sarcomeres in series and can undergo larger changes in length. This also means that a given change in muscle length of a longer muscle is accompanied by relatively smaller length changes at the sarcomere level, hence a smaller change in myofilament overlap and force produced. Because the dimensionless Force-length relationship (see Force-length) applies to behavior at the sarcomere level, it must be scaled by the length of the muscle to account for this architectural effect on the behavior of whole muscle. The morphometric parameter used for this purpose is optimal fascicle length, which is the distance between the aponeurosis sites measured along the orientation of the muscle fibers, measured at the length for which the muscle produces its maximal isometric force. Fascicle length is measured instead of muscle belly length because muscle fibers in many muscles are oriented obliquely (see pennation angle) with respect to the long axis of the muscle belly.

### Muscle mass ($$m$$)

The amount of force that a muscle can produce depends on the total number of myofilaments arranged in parallel. Because myofilaments are densely packed through muscle fibers, this number is proportional to the physiological cross-sectional area (PCSA) of the muscle fibers. Because those fibers may be arranged in a pennate form, the anatomical cross-section of the muscle may not be perpendicular to all the muscle fibers, so would not reflect their true cross-sectional area. More reliable measurements can be made of the mass of the muscle, which can be converted to its volume by dividing by density. This volume must be the product of the length of the contractile elements ($$L_{ce}$$), which can be measured reliably by dissecting fascicles from the muscle, assuming that they are of uniform length (common) and at optimal length $$L_{0}$$ (any discrepancy can be detected and corrected by observing sarcomere length under a microscope).

$PCSA = \frac{m}{\rho L_{ce}}$

Specific tension, ε, or maximal isometric force produced per unit of cross-sectional area, has been shown through single fiber (Lucas et al., 1987) and whole muscle (Spector et al., 1980, Brown, 1998) studies to be the same across different fiber types and was chosen to be 31.8 N/$$cm^2$$. Thus maximal isometric force ($$F_{0}$$) can be estimated using the following equation:

$F_0 = \varepsilon *PCSA$

Muscle fascicle length in the equation is set to $$L_{0}$$, because ε is defined at that length. The modeled muscle’s $$F_{0}$$ is used to scale the dimensionless relationships to obtain force output in absolute units.

### Fiber composition

Although maximal isometric force per unit area is roughly the same for different muscle fiber types, there are substantial differences in their contractile properties. Fiber types differ in terms of sarcomere length and velocity effects on force output, rise and fall dynamics of contraction, and especially energy consumption (see Energy related to cross-bridges). The recruitment order of motor units also depends on their fiber type (see Frequency-recruitment). It is therefore important that different fiber types are modeled independently and that the fiber composition of the modeled muscle is known accurately.

Determining fiber composition is often challenging because of the many classifications of fiber types and because of the difficulty in measuring the attributes that distinguish them. Fibers differ mainly in terms of their myosin isoforms (namely the myosin light and heavy chains) and the size of their sarcoplasmic reticulum and transverse tubule system that together influence force generation under dynamic conditions. There are also large differences in mitochondrial content and blood supply, which affect the metabolic cost of muscle contraction as well as the fiber’s susceptibility to fatigue. All of these tend to covary in healthy muscle fibers, but may become heterogenous in muscle fibers that are diseased or undergoing trophic changes as a result of abnormal or rapidly changing usage patterns. Because differences in myosin heavy chain isoform and mitochondrial content have the largest physiological effects, they are used most commonly to classify fiber types. Slow-twitch (type 1) fibers have a different distribution of myosin heavy chain isoforms than fast-twitch (type 2) fibers that leads to slower cross-bridge kinetics, as evidenced by the relatively slow rate at which they hydrolyze ATP. Type 1 fibers and some type 2 fibers (namely type 2a) have a high volume fraction of mitochondria, which can be assessed indirectly by measuring the concentrations of enzymes such as succinate dehydroginase, for example, which is normally present in mitochondria. These histochemical approaches require a tissue sample, which could be too invasive to obtain from humans. It is also worth noting that the results may be biased depending on the location in which it is obtained. Fast-twitch fibers tend to reside in the outer portion in muscle, so if a single sample is obtained from this location the measurement would underestimate the percentage of slow-twitch muscle. There are also noninvasive methods for measuring fiber composition such as magnetic resonance spectroscopy but to date they are relatively less accurate.

### Characteristic firing rate of muscle fiber type ($$f_{0.5}$$)

Muscles typically fire at rates ranging from 0.5$$f_{0.5}$$ to 2$$f_{0.5}$$, where $$f_{0.5}$$ is the firing rate at which muscle produces half of its maximal isometric force when its length is held at $$L_{0}$$ (see Optimal fascicle length). $$f_{0.5}$$ is an important parameter because cross-bridge activation relationships of different muscles that are normalized by this value become congruent and therefore a generic relationship can be obtained. This is useful because only this parameter is needed to construct the entire cross-bridge activation relationship for an arbitrary muscle to a high degree of accuracy. The actual firing rate behavior of motor units has been difficult to record, particularly at the higher levels of recruitment, so this useful normalization remains contentious.

### Tendon + aponeurosis length ($$L_{0}^T$$)

The tendon plus aponeurosis is composed of tightly packed strands of collagen that are oriented in parallel with the long axis of the musculotendon. As in muscle, the amount a tendon can stretch depends on its length and the amount of force generates in response to a stretch depends on its cross-sectional area (see Series elastic element). Longer tendons are composed of more collagen molecules arranged in series, each having similar material properties and dimensions, i.e. experiencing the same change in length in response to the same force. Longer tendons can therefore experience larger changes in length and will produce lower forces in response to the same stretch. The generic relationship of the series elastic element reflects the material properties of the constituent collagen fibers and should be scaled by the length of the tendon plus aponeurosis to be modeled. The length should be measured when muscle is exerting maximal isometric force on the tendon ($$L_{0}^T$$). If only the slack length of the tendon is known, then $$L_{0}^T$$ can be estimated as 105% of tendon slack length.

### Maximum musculotendon path length ($$L_{max}^{MT}$$)

As mentioned in Parallel elastic element, the amount of passive force generated by muscle at a particular length correlates with its maximum anatomical length. Maximum anatomical length is defined as the maximum length a muscle can have across all anatomical configurations of the joints it crosses. This can be estimated by first measuring maximum anatomical length of the musculotendon and subtracting the length corresponding to the tendon ($$L_{0}^T$$).

### Pennation angle

Muscle fibers of many muscles are oriented obliquely with respect to the tendon and so not all of the force that they generate is transmitted to the bone segments. It is reduced by the cosine of the pennation angle, with the remainder of the contractile force producing hydrostatic compression of the muscle. Even for relatively large pennation angles, however, this loss in total force transmission in isometric muscle is relatively low (Scott and Winter, 1991), so this model does not account for this effect and does not require specification of pennation angle. Modeling the effects of pennation angle would be important for the relatively uncommon scenario in which it is large and changes rapidly by large amounts during a contraction, as it would alter the kinematics of the muscle fascicles that have a strong influence on force production.

## References

Balnave CD, Allen DG (1996) The effect of muscle length on intracellular calcium and force in single fibres from mouse skeletal muscle. Journal of Physiology 492.3:705-713.

Barclay CJ (1996) Mechanical efficiency and fatigue of fast and slow muscles of the mouse. Journal of Physiology 497 (Pt 3):781-794.

Barclay CJ, Lichtwark GA, Curtin NA (2008) The energetic cost of activation in mouse fast-twitch muscle is the same whether measured using reduced filament overlap or N-benzyl-p-toluenesulphonamide. Acta Physiol 193:381-391.

Barclay CJ, Woledge RC, Curtin NA (2010a) Inferring crossbridge properties from skeletal muscle energetics. Prog Biophys Mol Biol 102:53-71.

Barclay CJ, Woledge RC, Curtin NA (2010b) Is the efficiency of mammalian (mouse) skeletal muscle temperature dependent? J Physiol 588:3819-3831.

Brown IE (1998) Measured and Modeled Properties of Mammalian Skeletal Muscle. vol. Ph.D: Queen's University.

Brown IE, Cheng EJ, Loeb GE (1999) Measured and modeled properties of mammalian skeletal muscle. II. The effects of stimulus frequency on force-length and force-velocity relationships. Journal of Muscle Research and Cell Motility 20:643.

Brown IE, Kim DH, Loeb GE (1998) The effect of sarcomere length on triad location in intact feline caudofeomoralis muscle fibres. JMuscle ResCell Motil 19:473-477.

Brown IE, Liinamaa TL, Loeb GE (1996a) Relationships between range of motion, L0, and passive force in five strap-like muscles of the feline hind limb. JMorphol 230:69-77.

Brown IE, Loeb GE (1998) Post-activation potentiation - A clue for simplifying models of muscle dynamics. Am Zool 38:743-754.

Brown IE, Loeb GE (1999a) Measured and modeled properties of mammalian skeletal muscle. I. The effects of post-activation potentiation on the time course and velocity dependencies of force production. JMuscle ResCell Motil 20:443-456.

Brown IE, Loeb GE (1999b) Measured and modeled properties of mammalian skeletal muscle: III. The effects of stimulus frequency on stretch-induced force enhancement and shortening-induced force depression. J Muscle Res Cell Motility.

Brown IE, Loeb GE (2000) Measured and modeled properties of mammalian skeletal muscle: IV. dynamics of activation and deactivation. JMuscle ResCell Motil 21:33-47.

Brown IE, Scott SH, Loeb GE (1996b) Mechanics of feline soleus: II. Design and validation of a mathematical model. JMuscle ResCell Motil 17:221-233.

Brunello E, Reconditi M, Elangovan R, Linari M, Sun Y, Narayanan T, Panine P, Piazzesi G, Irving M, Lombardi V (2007) Skeletal muscle resists stretch by rapid binding of the second motor domain of myosin to actin. PNAS 104:20114-20119.

Burke RE, Levine DN, Tsairis P, Zajac FE (1973) Physiological types and histochemical profiles in motor units of the cat gastrocnemius. Journal of Physiology 234:723-748.

Cannell M, Allen D (1984) Model of calcium movements during activation in the sarcomere of frog skeletal muscle. Biophysical Journal 45:913-925.

Colombini B, Nocella M, Benelli G, Cecchi G, Bagni MA (2008) Effect of temperature on cross-bridge properties in intact frog muscle fibers. Am J Physiol Cell Physiol 294:C1113-C1117.

Diamant J, Keller A, Baer E, Litt M, Arridge R (1972) Collagen; ultrastructure and its relation to mechanical properties as a function of ageing. Proceedings of the Royal Society of London Series B Biological Sciences 180:293-315.

Epstein M, Wong M, Herzog W (2006) Should tendon and aponeurosis be considered in series? J Biomech 39:2020-2025.

Ford LE, Huxley AF, Simmons RM (1977) Tension responses to sudden length change in stimulated frog muscle fibres near slack length. Journal of Physiology 269:441-515.

Ford LE, Huxley AF, Simmons RM (1981) The relation between stiffness and filament overlap in stimulated frog muscle fibres. Journal of Physiology 311:219-249.

Ford LE, Huxley AF, Simmons RM (1985) Tension transients during steady shortening of frog muscle fibres. Journal of Physiology 361:131-150.

Fusi L, Reconditi M, Linari M, Brunello E, Elangovan R, Lombardi V, Piazzesi G (2010) The mechanism of the resistance to stretch of isometrically contracting single muscle fibres. J Physiol 588:495-510.

Gasser HS, Hill AV (1924) The dynamics of muscular contraction. Proc R Soc Lond 96:398-437.

Gordon AM, Homsher E, Regnier M (2000) Regulation of contraction in striated muscle. Physiol Rev 80:853-924.

Gordon AM, Huxley AF, Julian FJ (1966) The variation in isometric tension with sarcomere length in vertebral muscle fibres. Journal of Physiology 184:170-192.

Henneman E, Mendell LM (1981) Functional organization of the motoneuron pool and its inputs. In: Handbook of Physiology Sect I The Nervous System Vol II, Part 1(Brooks, V. B., ed), pp 423-507 Washington, DC: American Physiological Society.

Hill AV (1938) The heat of shortening and the dynamic constants of muscle. Proceedings of the Royal Society of London 126:136-195.

Hoffer JA, Caputi AA, Pose IE, Griffiths RI (1989) Roles of muscle activity and load on the relationship between muscle spindle length and whole muscle length in the freely walking cat. ProgBrain Res 80:75-85.

Hughes JR (1958) Post-tetanic potentiation. Physiological Reviews 38:91-113.

Huxley AF (1957) Muscle structure and theories of contraction. Progress in Biophysics and Molecular Biology 7:255-318.

Irving T, Wu Y, Bekyarova T, Farman GP, Fukuda N, Granzier H (2011) Thick-filament strain and interfilament spacing in passive muscle: effect of titin-based passive tension. Biophys J 100:1499-1508.

Joyce GS, Rack PMH, Westbury DR (1969) Mechanical properties of cat soleus muscle during controlled lengthening and shortening movements. J Physiol (Lond) 204:461-474.

Kawakami Y, Lieber RL (2000) Interaction between series compliance and sarcomere kinetics determines internal sarcomere shortening during fixed-end contraction. J Biomech 33:1249-1255.

Kaya M, Higuchi H (2010) Nonlinear elasticity and an 8-nm working stroke of single myosin molecules in myofilaments. Science 329:686-689.

Leijendekker WJ, Elzinga G (1990) Metabolic recovery of mouse extensor digitorum longus and soleus muscle. Pflügers Archiv European Journal of Physiology 416:22-27.

Loeb GE, Brown IE, Lan N, Davoodi R (2002) The importance of biomechanics. AdvExpMedBiol 508:481-487.

Loeb GE, Richmond FJR (1989) Motor partitioning: epiphenomena masquerading as control theory. Behav Brain Sci 12:660-661.

Lucas SM, Ruff RL, Binder MD (1987) Specific tension measurements in single soleus and medial gastrocnemius muscle fibres of the cat. Experimental Neurology 95:142-154.

Ma S, Zahalak GI (1991) A distribution-moment model of energetics in skeletal muscle. J Biomech 24:21-35.

Magid A, Law DJ (1985) Myofibrils bear most of the resting tension in frog skeletal muscle. Science 230:1280-1282.

Mansson A (2010) Significant impact on muscle mechanics of small nonlinearities in myofilament elasticity. Biophys J 99:1869-1875.

Matthews PB (1981) Muscle spindles: their messages and their fusimotor supply. Comprehensive Physiology.

Monster AW, Chan H (1977) Isometric force production by motor units of extensor digitorum communis muscle in man. Journal of Neurophysiology 40:1432-1443.

Nocella M, Bagni MA, Cecchi G, Colombini B (2013) Mechanism of force enhancement during stretching of skeletal muscle fibres investigated by high time-resolved stiffness measurements. J Muscle Res Cell Motil 34:71-91.

Perreault EJ, Heckman CJ, Sandercock TG (2003) Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. Journal of Biomechanics 36:211-218.

Piazzesi G, Francini F, Linari M, Lombardi V (1992) Tension transients during steady lengthening of tetanized muscle fibres of the frog. J Physiol 445:659-711.

Piazzesi G, Reconditi M, Linari M, Lucii L, Bianco P, Brunello E, Decostre V, Stewart A, Gore DB, Irving TC, Irving M, Lombardi V (2007) Skeletal Muscle Performance Determined by Modulation of Number of Myosin Motors Rather Than Motor Force or Stroke Size. Cell 131:784-795.

Rack PMH, Westbury DR (1969) The effects of length and stimulus rate on tension in the isometric cat soleus muscle. Journal of Physiology 204:443-460.

Rall JA, Woledge RC (1990) Influence of temperature on mechanics and energetics of muscle contraction Am J Physiol Regul Integr Comp Physiol 259:197-203.

Rassier DE, Macintosh BR (2002) Length-dependent twitch contractile characteristics of skeletal muscle. Can J Physiol Pharmacol 80:993-1000.

Scott SH, Brown IE, Loeb GE (1996) Mechanics of feline soleus: I. Effect of fascicle length and velocity on force output. JMuscle ResCell Motil 17:205-218.

Scott SH, Loeb GE (1995) Mechanical properties of aponeurosis and tendon of the cat soleus muscle during whole-muscle isometric contractions. JMorphol 224:73-86.

Scott SH, Winter DA (1991) A comparison of three muscle pennation assumptions and their effect on isometric and isotonic force. J Biomech 24:163-167.

Smith IC, Gittings W, Huang J, McMillan EM, Quadrilatero J, Tupling AR, Vandenboom R (2013) Potentiation in mouse lumbrical muscle without myosin light chain phosphorylation: Is resting calcium responsible? The Journal of general physiology 141:297-308.

Spector SA, Gardiner PF, Zernicke RF, Roy RR, Edgerton VR (1980) Muscle architecture and the force velocity characteristics of the cat soleus and medial gastrocnemius: Implications for motor control. Journal of Neurophysiology 44:951-960.

Stull JT, Kamm KE, Vandenboom R (2011) Myosin light chain kinase and the role of myosin light chain phosphorylation in skeletal muscle. Arch Biochem Biophys 510:120-128.

Thelen DG (2003) Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults. Journal of Biomechanical Engineering-Transactions of the Asme 125:70-77.

Tilp M, Steib S, Herzog W (2012) Length changes of human tibialis anterior central aponeurosis during passive movements and isometric, concentric, and eccentric contractions. Eur J Appl Physiol 112:1485-1494.

Tsianos GA, Rustin C, Loeb GE (2012) Mammalian muscle model for predicting force and energetics during physiological behaviors. IEEE Trans Neural Syst Rehabil Eng 20:117-133.

Wang K, McCarter R, Wright J, Beverly J, Ramirez-Mitchell R (1993) Viscoelasticity of the sarcomere matrix of skeletal muscles. The titin-myosin composite filament is a dual-stage molecular spring. Biophysical Journal 64:1161-1177.

Zahalak GI (1981) A Distribution-Moment Approximation for Kinetic Theories of Muscular Contraction. Mathematical Biosciences 55:89-114.

Zajac FE (1989) Muscle and tendon: properties, models, scaling and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering 17:359-411.

Zajac FE, Zomlefer MR, Levine WS (1981) Hindlimb muscular activity, kinetics and kinematics of cats jumping to their maximum achievable heights. J Exp Biol 91:73-86.