Abstract
The force developed by actively lengthened muscle depends on different structures across different scales of lengthening. For small perturbations, the active response of muscle is well captured by a linear-time-invariant (LTI) system: a stiff spring in parallel with a light damper. The force response of muscle to longer stretches is better represented by a compliant spring that can fix its end when activated. Experimental work has shown that the stiffness and damping (impedance) of muscle in response to small perturbations is of fundamental importance to motor learning and mechanical stability, while the huge forces developed during long active stretches are critical for simulating and predicting injury. Outside of motor learning and injury, muscle is actively lengthened as a part of nearly all terrestrial locomotion. Despite the functional importance of impedance and active lengthening, no single muscle model has all of these mechanical properties. In this work, we present the viscoelastic-crossbridge active-titin (VEXAT) model that can replicate the response of muscle to length changes great and small. To evaluate the VEXAT model, we compare its response to biological muscle by simulating experiments that measure the impedance of muscle, and the forces developed during long active stretches. In addition, we have also compared the responses of the VEXAT model to a popular Hill-type muscle model. The VEXAT model more accurately captures the impedance of biological muscle and its responses to long active stretches than a Hill-type model and can still reproduce the force-velocity and force-length relations of muscle. While the comparison between the VEXAT model and biological muscle is favorable, there are some phenomena that can be improved: the low frequency phase response of the model, and a mechanism to support passive force enhancement.
eLife assessment
This is a valuable study that develops a new model of the way muscle responds to perturbations, synthesizing models of how it responds to small and large perturbations, both of which are used to predict how muscles function for stability but also how they can be injured, and which tend to be predicted poorly by classic Hill-type models. The evidence presented to support the model is solid, since it outperforms Hill-type models in a variety of conditions. Although the combination of phenomenological and mechanistic aspects of the model may sometimes make it challenging to interpret the output, the work will be of interest to those developing realistic models of the stability and control of movement in humans or other animals.
1 Introduction
The stiffness and damping of muscle are properties of fundamental importance for motor control, and the accurate simulation of muscle force. The central nervous system (CNS) exploits the activation-dependent stiffness and damping (impedance) of muscle when learning new movements [1], and when moving in unstable [2] or noisy environments [3]. Reaching experiments using haptic manipulanda show that the CNS uses co-contraction to increase the stiffness of the arm when perturbed by an unstable force field [4]. With time and repetition, the force field becomes learned and co-contraction is reduced [1].
The force response of muscle is not uniform, but varies with both the length and time of perturbation. Under constant activation and at a consistent nominal length, Kirsch et al. [5] were able to show that muscle behaves like a linear-time-invariant (LTI) system in response to small1 perturbations: a spring-damper of best fit captured over 90% of the observed variation in muscle force for small perturbations (1-3.8% optimal length) over a wide range of bandwidths (4-90Hz). When active muscle is stretched appreciably, titin can develop enormous forces [7], [8], which may prevent further lengthening and injury. The stiffness that best captures the response of muscle to the small perturbations of Kirsch et al. [5] is far greater than the stiffness that best captures the response of muscle to large perturbations [7], [8]. Since everyday movements are often accompanied by both large and small kinematic perturbations, it is important to accurately capture these two processes.
However, there is likely no single muscle model that can replicate the force response of muscle to small [5] and large perturbations [7], [8] while also retaining the capability to reproduce the experiments of Hill [9] and Gordon et al. [10]. Unfortunately, this means that simulation studies that depend on an accurate representation of muscle impedance may reach conclusions well justified in simulation but not in reality. In this work, we focus on formulating a mechanistic muscle model2 that can replicate the force response of active muscle to length perturbations both great and small.
There are predominantly three classes of models that are used to simulate musculoskeletal responses: phenomenological models constructed using Hill’s famous force-velocity relationship [9], mechanistic Huxley [11]–[13] models in which individual elastic crossbridges are incorporated, and linearized muscle models [14], [15] which are accurate for small changes in muscle length. Kirsch et al. [5] demon-strated that, for small perturbations, the force response of muscle is well represented by a spring in parallel with a damper. Neither Hill nor Huxley models are likely to replicate Kirsch et al.’s [5] experiments because a Hill muscle model [16], [17] does not contain any active spring elements; while a Huxley model lacks an active damping element. Although linearized muscle models can replicate Kirsch et al.’s experiment [5], these models are only accurate for small changes in length and cannot replicate the Hill’s nonlinear force-velocity relation [9], nor Gordon et al.’s [10] nonlinear force-length relation. However, there have been significant improvements to the canonical forms of phenomenological, mechanistic, and linearized muscle models that warrant closer inspection.
Several novel muscle models have been proposed to improve upon the accuracy of Hill-type muscle models during large active stretches. Forcinito et al. [18] modeled the velocity dependence of muscle using a rheological element3 and an elastic rack rather than embedding the force-velocity relationship in equations directly, as is done in a typical Hill model [16], [17]. This modification allows Forcinito et al.’s [18] model to more faithfully replicate the force development of active muscle, as compared to a Hill-type model, during ramp length changes of ≈ 10%4 of the optimal CE length, and across velocities of 4 − 11% of the maximum contraction velocity5. Haeufle et al. [20] made use of a serial-parallel network of spring-dampers to allow their model to reproduce Hill’s force-velocity relationship [9] mechanistically rather than embedding the experimental curve directly in their model. This modification allowed Haeufle et al.’s model to simulate high speed reaching movements that agree more closely with experimental data [21] than is possible with a typical Hill model. Gü nther et al. [22] evaluated how accurately a variety of spring-damper models were able to reproduce the microscopic increases in crossbridge force in response to small length changes. While each of these models improves upon the force response of the Hill model to ramp length changes, none are likely to reproduce Kirsch et al.’s experiment [5] because the linearized versions of these models lead to a serial, rather than a parallel, connection of a spring and a damper: Kirsch et al. [5] specifically showed (see Figure 3 of [5]) that a serial connection of a spring-damper fails to reproduce the phase shift between force and length present in their experimental data.
Titin [23], [24] has been more recently investigated to explain how lengthened muscle can develop active force when lengthened both within, and beyond, actin-myosin overlap [8]. Titin is a gigantic multi-segmented protein that spans a half-sarcomere, attaching to the Z-line at one end and the middle of the thick filament at the other end [25]. In skeletal muscle, the two sections nearest to the Z-line, the proximal immunoglobulin (IgP) segment and the PEVK segment — rich in the amino acids proline (P), glutamate (E), valine (V) and lysine (K) — are the most compliant [26] since the distal immunoglobulin (IgD) segments bind strongly to the thick filament [27]. Titin has proven to be a complex filament, varying in composition and geometry between different muscle types [28], [29], widely between species [30], and can apply activation dependent forces to actin [31]. It has proven challenging to determine which interactions dominate between the various segments of titin and the other filaments in a sarcomere. Experimental observations have reported titin-actin interactions at myosin-actin binding sites [32], [33], between titin’s PEVK region and actin [34], [35], between titin’s N2A region and actin [36], and between the PEVK-IgD regions of titin and myosin [37]. This large variety of experimental observations has led to a correspondingly large number of proposed hypotheses and models, most of which involve titin interacting with actin [38]–[43], and more recently with myosin [44].
The addition of a titin element to a model will result in more accurate force production during large active length changes, but does not affect the stiffness and damping of muscle at modest sarcomere lengths because of titin’s relatively low stiffness. At sarcomere lengths of
Although most motor control simulations [2], [45]–[48] make use of the canonical linearized muscle model, phenomenological muscle models have also been used and modified to include stiffness. Sartori et al. [49] modeled muscle stiffness by evaluating the partial derivative of the force developed by a Hill-type muscle model with respect to the contractile element (CE) length. Although this approach is mathematically correct, the resulting stiffness is heavily influenced by the shape of the force-length curve and can lead to inaccurate results: at the optimal CE length this approach would predict an active muscle stiffness of zero since the slope of the force-length curve is zero; on the descending limb this approach would predict a negative active muscle stiffness since the slope of the force-length curve is negative. In contrast, CE stiffness is large and positive near the optimal length [5], and there is no evidence for negative stiffness on the descending limb of the force-length curve [7]. Although the stiffness of the CE can be kept positive by shifting the passive force-length curve, which is at times used in finite-element-models of muscle [50], this introduces a new problem: the resulting passive CE stiffness cannot be lowered to match a more flexible muscle. In contrast, De Groote et al. [51], [52] modeled short-range-stiffness using a stiff spring in parallel with the active force element of a Hill-type muscle model. While the approach of De Groote et al. [51], [52] likely does improve the response of a Hill-type muscle model for small perturbations, there are several drawbacks: the short-range-stiffness of the muscle sharply goes to zero outside of the specified range whereas in reality the stiffness is only reduced [5] (see Fig. 9A); the damping of the canonical Hill-model has been left unchanged and likely differs substantially from biological muscle [5].
In this work, we propose a model that can capture the force development of muscle to perturbations that vary in size and timescale, and yet is described using only a few states making it well suited for large-scale simulations. When active, the response of the model to perturbations within actin-myosin overlap is dominated by a viscoelastic crossbridge element that has different dynamics across time-scales: over brief time-scales the viscoelasticity of the lumped crossbridge dominates the response of the muscle [5], while over longer time-scales the force-velocity [9] and force-length [10] properties of muscle dominate. To capture the active forces developed by muscle beyond actin-myosin overlap we added an active titin element which, similar to existing models [38], [40], features an activation-dependent6 interaction between titin and actin. To ensure that the various parts of the model are bounded by reality, we have estimated the physical properties of the viscoelastic crossbridge element as well as the active titin element using data from the literature.
While our main focus is to develop a more accurate muscle model, we would like the model to be well suited to simulating systems that contain tens to hundreds of muscles. Although Huxley models have been used to simulate whole-body movements such as jumping [53], the memory and processing requirements associated with simulating a single muscle with thousands of states is high. Instead of modeling the force development of individual crossbridges, we lump all of the crossbridges in a muscle together so that we have a small number of states to simulate per muscle.
To evaluate the proposed model, we compare simulations of experiments to original data. We examine the response of active muscle to small perturbations over a wide band-width by simulating the stochastic perturbation experiments of Kirsch et al. [5]. Herzog et al.’s [7] active-lengthening experiments are used to evaluate the response of the model when it is actively lengthened within actin-myosin overlap. Next, we use Leonard et al.’s [8] active lengthening experiments to see how the model compares to reality when it is actively lengthened beyond actin-myosin overlap. In addition, we examine how well the model can reproduce the force-velocity experiments of Hill [9] and force-length experiments of Gordon et al. [10]. Since Hill-type models are so commonly used, we also replicate all of the simulated experiments using Millard et al.’s [17] Hill-type muscle model to make the differences between these two types of models clear.
2 Model
We begin by treating whole muscle as a scaled half-sarcomere that is pennated at an angle α with respect to a tendon (Fig. 1A). The assumption that mechanical properties scale with size is commonly used when modeling muscle [16] and makes it possible to model vastly different musculotendon units (MTUs) by simply changing the architectural and contraction properties: the maximum isometric force
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig1.tif/full/max/0/default.jpg)
The name of the VEXAT model comes from the viscoelastic crossbridge and active titin elements (A.) in the model. Active tension generated by the lumped crossbridge flows through actin, myosin, and the adjacent sarcomeres to the attached tendon (B.). Titin is modeled as two springs of length 𝓁 1 and 𝓁 2 in series with the rigid segments LT12 and LIgD. Viscous forces act between titin and actin in proportion to the activation of the muscle (C.), which reduces to negligible values in a purely passive muscle (D.). We modeled actin and myosin as rigid elements; the XE, titin, and the tendon as viscoelastic elements; and the ECM as an elastic element.
The proposed model has several additional properties that we assume scale with
where
where
To reduce the number of states needed to simulate the VEXAT model, we lump all of the attached crossbridges into a single lumped crossbridge element (XE) that attaches at 𝓁 S (Fig. 1A) and has intrinsic stiffness and damping properties that vary with the activation and force-length properties of muscle. The active force developed by the XE at the attachment point to actin is transmitted to the main myosin filament, the M-line, and ultimately to the tendon (Fig. 1B). In addition, since the stiffness of actin [60] and myosin filaments [61] greatly exceeds that of crossbridges [62], we treat actin and myosin filaments as rigid to reduce the number of states needed to simulate this model. Similarly, we have lumped the six titin filaments per half-sarcomere (Fig. 1A) together to further reduce the number of states needed to simulate this model.
The addition of a titin filament to the model introduces an additional active load-path (Fig. 1C) and an additional passive load-path (Fig. 1D). As is typical [16], [17], we assume that the passive elasticity of these structures scale linearly with
As previously mentioned, there are several theories to explain how titin interacts with the other filaments in activated muscle. While there is evidence for titin-actin interaction near titin’s N2A region [36], there is also support for a titin-actin interaction occurring near titin’s PEVK region [34], [35], and for a titin-myosin interaction near the PEVK-IgD region [37]. For the purposes of our model, we will assume a titin-actin interaction because current evidence weighs more heavily towards a titin-actin interaction than a titin-myosin interaction. Next, we assume that the titin-actin interaction takes place somewhere in the PEVK segment for two reasons: first, there is evidence for a titin-actin interaction [34], [35] in the PEVK segment; and second, there is evidence supporting an interaction at the proximal end of the PEVK segment (N2A-actin interaction) [36]. We have left the point within the PEVK segment that attaches to actin as a free variable since there is some uncertainty about what part of the PEVK segment interacts with actin. The nature of the mechanical interaction between titin and the other filaments in an active sarcomere remains uncertain. Here we assume that this interaction is not a rigid attachment, but instead is an activation dependent damping to be consistent with the observations of Kellermayer and Granzier [31] and Dutta et al. [36]: adding titin filaments and calcium slowed, but did not stop, the progression of actin filaments across a plate covered in active crossbridges (an in-vitro motility assay). When activated, we assume that the amount of damping between titin and actin scales linearly with
After lumping all of the crossbridges and titin filaments together we are left with a rigid-tendon MTUmodel that has two generalized positions
and an elastic-tendon MTUmodel that has three generalized positions
Given these generalized positions, the path length 𝓁 P, and a pennation model, all other lengths in the model can be calculated. Here we use a constant thickness
pennation model to evaluate the pennation angle
of a CE with a rigid-tendon, and
to evaluate the pennation angle of a CE with an elastic-tendon. We have added a small compressive element KE (Fig. 1A) to prevent the model from reaching the numerical singularity that exists as
of an elastic-tendon model is the difference between the path length and the CE length along the tendon. The length of the XE
is the difference between the half-sarcomere length and the sum of the average point of attachment 𝓁 S and the length of the myosin filament LM. The length of 𝓁 2, the lumped PEVK-IgD segment, is
the difference between the half-sarcomere length and the sum of the length from the Z-line to the actin binding site on titin (𝓁 1) and the length of the IgD segment that is bound to myosin (LIgD). Finally, the length of the extra-cellular-matrix 𝓁 ECM is simply
half the length of the CE since we are modeling a half-sarcomere.
We have some freedom to choose the state vector of the model and the differential equations that define how the muscle responds to length and activation changes. The experiments we hope to replicate depend on phenomena that take place at different time-scales: Kirsch et al.’s [5] stochastic perturbations evolve over brief time-scales, while all of the other experiments take place at much longer time-scales. Here we mathematically decouple phenomena that affect brief and long time-scales by making a second-order model that has states of the average point of crossbridge attachment 𝓁 S, and velocity v S. When the activation a state and the titin-actin interaction model are included, the resulting rigid-tendon model that has a total of four states
and the elastic-tendon model has
five states. For the purpose of comparison, a Hill-type muscle model with a rigid-tendon has a single state (a), while an elastic-tendon model has two states (a and 𝓁 M) [17].
Before proceeding, a small note on notation: throughout this work we will use an underbar to indicate a vector, bold font to indicate a curve, a tilde for a normalized quantity, and a capital letter to indicate a constant. Unless indicated otherwise, curves are constructed using C2 continuous7 Bézier splines so that the model is compatible with gradient-based optimization. Normalized quantities within the CE follow a specific convention: lengths and velocities are normalized by the optimal CE length
To evaluate the state derivative of the model, we require equations for
over longer time-scales, where f L(·) is the active-force-length curve (Fig. 2A), f PE(·) is the passive-force-length curve (Fig. 2A), and f V(·) is the force-velocity (Fig. 2B).
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig2.tif/full/max/0/default.jpg)
The model relies on Bézier curves to model the nonlinear effects of the active-force-length curve, the passive-force-length curves (A.), and the force-velocity curve (B.). Since nearly all of the reference experiments used in Sec. 3 have used cat soleus, we have fit the active-force-length curve (f L(·)) and passive-force-length curves (f PE(·)) to the cat soleus data of Herzog and Leonard 2002 [7]. The concentric side of the force-velocity curve (f V(·)) has been fitted to the cat soleus data of Herzog and Leonard 1997 [63].
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig3.tif/full/max/0/default.jpg)
The passive force-length curve has been decomposed such that 56% of it comes from the ECM while 44% comes from titin to match the average of ECM-titin passive force distribution (which ranges from 43%-76%) reported by Prado et al. [58] (A.). The elasticity of the titin segment has been further decomposed into two serially connected sections: the proximal section consisting of the T12, proximal IgP segment and part of the PEVK segment, and the distal section consisting of the remaining PEVK section and the distal Ig segment (B.). The stiffness of the IgP and PEVK segments has been chosen so that the model can reproduce the movements of IgP/PEVK and PEVK/IgD boundaries that Trombitás et al. [26] (C.) observed in their experiments. The curves that appear in subplots A. and B. come from scaling the two-segmented human soleus titin model to cat soleus muscle. The curves that appear in subplot C compare the human soleus titin model’s IgP, PEVK, and IgD force-length relations to the data of Trombitás et al. [26] (see in Appendix B for details).
The normalized tension developed by the VEXAT model
differs from that of a Hill model, Eqn. 14, because it has no explicit dependency on
We set the first term of
where τ S is a time constant and
The remaining two terms,
The passive force developed by the CE in Eqn. 15 is the sum of the elastic forces (Fig. 3A) developed by the force-length curves of titin
The curves that form
When active muscle is lengthened, it produces an enhanced force that persists long after the lengthening has ceased called residual force enhancement (RFE) [7]. For the purposes of the VEXAT model, we assume that RFE is produced by titin. Experiments have demonstrated RFE on both the ascending limb [66] and descending limb of the force-length [7] relation. The amount of RFE depends both on the final length of the stretch [67] and the magnitude of the stretch: on the ascending limb the amount of RFE varies with the final length but not with stretch magnitude, while on the descending limb RFE varies with stretch magnitude. To develop RFE, we assume that some point of the PEVK segment bonds with actin through an activation-dependent damper. The VEXAT model’s distal segment of titin, 𝓁 2, can contribute to RFE when the titin-actin bond is formed and CE is lengthened beyond
that transitions from zero to one as
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl1.tif/full/max/0/default.jpg)
The VEXAT and Hill model’s elastic-tendon cat soleus MTU parameters. The VEXAT model uses all of the Hill model’s parameters which are highlighted in grey. Short forms are used to indicate: length ‘len’, velocity ‘vel’, acceleration ‘acc’, half ‘h’, activation ‘act’, segment ‘seg’, threshold ‘thr’, and stiffness ‘stiff’. The letters ‘R’ or ‘H’ in front of a reference mean the data is from a rabbit or a human, otherwise the data is from cat soleus. The letters following a reference indicate how the data was used to create the parameter: ‘C’ calculated, ‘F’ fit, ‘E’ estimated, or ‘S’ scaled. Most of the VEXAT model’s XE and titin parameters can be used as rough parameter guesses for other muscles because we have expressed these parameters in a normalized space: the values will scale appropriately with changes to
At very long CE lengths, the modeled titin-actin bond can literally slip off of the end of the actin filament (Fig. 1A) when the distance between the Z-line and the bond,
The step-down function uL transitions from one to zero when the titin-actin bond
The strength of the titin-actin bond also appears to vary nonlinearly with activation. Fukutani and Herzog [68] observed that the absolute RFE magnitude produced by actively lengthened fibers is similar between normal and reduced contractile force states. Since these experiments [68] were performed beyond the optimal CE length, titin could be contributing to the observed RFE as previously described. The consistent pattern of absolute RFE values observed by Fukutani and Herzog [68] could be produced if the titin-actin bond saturated at its maximum strength even at a reduced contractile force state. To saturate the titin-actin bond, we use a final smooth step function
where A° is the threshold activation level at which the bond saturates. While we model the strength of the titin-actin bond as being a function of activation, which is proportional Ca2+ concentration [69], this is a mathematical convenience. The work of Leonard et al. [8] makes it clear that both Ca2+ and crossbridge cycling are needed to allow titin to develop enhanced forces during active lengthening: no enhanced forces are observed in the presence of Ca2+ when crossbridge cycling is chemically inhibited. Putting this all together, the active damping acting between the titin and actin filaments is given by the product of
With a model of the titin-actin bond derived, we can focus on how the bond location moves in response to applied forces. Since we are ignoring the mass of the titin filament, the PEVK-attachment point is balanced by the forces applied to it and the viscous forces developed between titin and actin
due to the active
The assumption of whether the tendon is rigid or elastic affects how the state derivative is evaluated and how expensive it is to compute. While all of the position dependent quantities can be evaluated using Eqns. 6-11 and the generalized positions, evaluating the generalized velocities of a rigid-tendon and elastic-tendon model differ substantially.
The CE velocity v M and pennation angular velocity
Evaluating the CE rate of lengthening, v M, for an elastictendon muscle model is more involved. As is typical of lumped parameter muscle models [16], [17], [70], here we assume that difference in tension,
is negligible10. During our preliminary simulations it became clear that treating the tendon as an idealized spring degraded the ability of the model to replicate the experiment of Kirsch et al. [5] particularly at high frequencies. Kirsch et al. [5] observed a linear increase in the gain and phase profile between the output force and the input perturbation applied to the muscle. This pattern in gain and phase shift can be accurately reproduced by a spring in parallel with a damper. Due to the way that impedance combines in series 11, the models of both the CE and the tendon need to have parallel spring and damper elements so that the entire MTU, when linearized, appears to be a spring in parallel with a damping element. We model tendon force using a nonlinear spring and damper model
where the damping coefficient
allowing us to evaluate the final state derivative in
3 Biological Benchmark Simulations
In order to evaluate the model, we have selected three experiments that capture the responses of active muscle to small, medium, and large length changes. The small
The VEXAT model requires the architectural muscle parameters (
3.1 Stochastic Length Perturbation Experiments
In Kirsch et al.’s [5] in-situ experiment, the force response of a cat’s soleus muscle under constant stimulation was measured as its length was changed by small amounts. Kirsch et al. [5] applied stochastic length perturbations (Fig. 4A) to elicit force responses from the muscle (in this case a spring-damper Fig. 4B) across a broad range of frequencies (4-90 Hz) and across a range of small length perturbations
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig4.tif/full/max/0/default.jpg)
Evaluating a system’s gain and phase response begins by applying a pseudo-random input signal to the system and measuring its output (A). Both the input and output signals (A) are transformed into the frequency domain by expressing these signals as an equivalent sum of scaled and shifted sinusoids (simple example shown in B and C). Each individual input sinusoid is compared with the output sinusoid of the same frequency to evaluate how the system scales and shifts the input to the output (D). This process is repeated across all sinusoid pairs to produce a function that describes how an input sinusoid is scaled (E) and shifted (F) to an output sinusoid using only the measured data (A).
As long as muscle can be considered to be linear (a sinusoidal change in length produces a sinusoidal change in force), then system identification methods [78], [79] can be applied to extract a relationship between length x(t) and force y(t). We will give a brief overview of system identification methods here to make methods and results clearer. First, the time-domain signals (x(t) and y(t)) are transformed into an equivalent representation in the frequency-domain (X(s) and Y (s)) as a sum of scaled and shifted sine curves (Fig. 4B and C) using a Fourier transform [78].
In the frequency domain, we identify an LTI system of best fit H(s) that describes how muscle transforms changes in length into changes in force such that Y (s) = H(s) X(s). Next, we evaluate how H(s) scales the magnitude (gain) and shifts the timing (phase) of a sinusoid in X(s) into a sinusoid of the same frequency in Y (s) (Fig. 4D). This process is repeated across all frequency-matched pairs of input and output sinusoids to build a function of how muscle scales (Fig. 4E) and shifts (Fig. 4F) input length sinusoids into output force sinusoids. The resulting transformation turns two complicated time-domain signals (Fig. 4A) into a clear relationship in the frequency-domain that describes how muscle transforms length changes into force changes: a very slow (0Hz) length change will result in an output force that is scaled by 4.5 and is in phase (Fig. 4E and F), a 35 Hz sinusoidal length change will produce an output force that is scaled by 4.9 and leads the input signal by 24° (Fig. 4E and F), and frequencies between 0 Hz and 35 Hz will be between these two signals in terms of scaling and phase. These patterns of gain and phase can be used to identify a network of spring-dampers that is equivalent to the underlying linear system (the system in Fig. 4 A, E, and F is a 4.46 N/mm spring in parallel with a 0.0089 Ns/mm damper). Since experimental measurements often contain noise and small nonlinearities, the mathematical procedure used to estimate H(s) and the corresponding gain and phase profiles is more elaborate than we have described (see Appendix D for details).
Kirsch et al. [5] used system identification methods to identify LTI mechanical systems that best describes how muscle transforms input length waveforms to output force waveforms. The resulting LTI system, however, is only accurate when the relationship between input and output is approximately linear. Kirsch et al. [5] used the coherence squared, (Cxy)2, between the input and output to evaluate the degree of linearity: Y(s) is a linear transformation of X(s) at frequencies in which (Cxy)2 is near one, while Y (s) cannot be described as a linear function of X(s) at frequencies in which (Cxy)2 approaches zero. By calculating (Cxy)2 between the length perturbation and force waveforms, Kirsch et al. [5] identified the bandwidth in which the muscle’s response is approximately linear. Kirsch et al. [5] set the lower frequency of this band to 4 Hz, and Fig. 3 of Kirsch et al. [5] suggests that this corresponds to (Cxy)2 ≥ 0.67 though the threshold for (Cxy)2 is not reported. The upper frequency of this band was set to the cutoff frequency of the low-pass filter applied to the input (15, 35, or 90 Hz). Within this bandwidth, Kirsch et al. [5] compared the response of the specimen to several candidate models and found that a parallel spring-damper fit the muscle’s frequency response best. Next, they evaluated the stiffness and damping coefficients that best fit the muscle’s frequency response [5]. Finally, Kirsch et al. evaluated how much of the muscle’s time-domain response was captured by the spring-damper of best fit by evaluating the variance-accounted-for (VAF) between the two time-domain signals
Astonishingly, Kirsch et al. [5] found that a spring-damper of best fit has a VAF of between 78-99%12 when compared to the experimentally measured forces f EXP. By repeating this experiment over a variety of stimulation levels (using both electrical stimulation and the crossed-extension reflex) Kirsch et al. [5] showed that these stiffness and damping coefficients vary linearly with the active force developed by the muscle. Further, Kirsch et al. [5] repeated the experiment using perturbations that had a variety of lengths (0.4 mm, 0.8mm, and 1.6mm) and bandwidths (15Hz, 35Hz, and 90Hz) and observed a peculiar quality of muscle: the damping coefficient of best fit increases as the bandwidth of the perturbation decreases (See Figures 3 and 10 of Kirsch et al. [5] for details). Here we simulate Kirsch et al.’s experiment [5] to determine, first, the VAF of the VEXAT model and the Hill model in comparison to a spring-damper of best fit; second, to compare the gain and phase response of the models to biological muscle; and finally, to see if the spring-damper coefficients of best fit for both models increase with active force in a manner that is similar to the cat soleus that Kirsch et al. studied [5].
To simulate the experiments of Kirsch et al. [5] we begin by creating the 9 stochastic perturbation waveforms used in the experiment that vary in perturbation amplitude (0.4mm, 0.8mm, and 1.6mm) and bandwidth (0-15 Hz, 0-35 Hz, and 0-90 Hz)13. The waveform is created using a vector that is composed of random numbers with a range of [−1, 1] that begins and ends with a series of zero-valued padding points.
Next, a forward pass of a 2nd order Butterworth filter is applied to the waveform and finally the signal is scaled to the appropriate amplitude (Fig. 5). The muscle model is then activated until it develops a constant tension at a length of
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig5.tif/full/max/0/default.jpg)
The perturbation waveforms are constructed by generating a series of pseudo-random numbers, padding the ends with zeros, by filtering the signal using a 2nd order low-pass filter (wave forms with -3dB cut-off frequencies of 90 Hz, 35 Hz and 15 Hz appear in A.) and finally by scaling the range to the desired limit (1.6mm in A.). Although the power spectrum of the resulting signals is highly variable, the filter ensures that the frequencies beyond the -3dB point have less than half their original power (B.).
When coupled with an elastic-tendon, the 15 Hz perturbations show that neither model can match the VAF of Kirsch et al.’s analysis [5] (compare Fig. 6A to G), while at 90Hz the VEXAT model reaches a VAF of 89% (Fig. 6D) which is within the range of 78-99% reported by Kirsch et al. [5]. In contrast, the Hill model’s VAF at 90 Hz remains low at 58% (Fig. 6J). While the VEXAT model has a gain profile in the frequency-domain that closer to Kirsch et al.’s data [5] than the Hill model (compare Fig. 6B to H and E to K), both models have a greater phase shift than Kirsch et al.’s data [5] at low frequencies (compare Fig. 6C to I and F to L). The phase response of the VEXAT model to the 90 Hz perturbation (Fig. 6F) shows the consequences of Eqn. 16: at low frequencies the phase response of the VEXAT model is similar to that of the Hill model, while at higher frequencies the model’s response becomes similar to a spring-damper. This frequency dependent response is a consequence of the first term in Eqn. 16: the value of τ S causes the response of the model to be similar to a Hill model at lower frequencies and mimic a spring-damper at higher frequencies. Both models show the same perturbation-dependent phase-response, as the damping coefficient of best fit increases as the perturbation bandwidth decreases: compare the damping coefficient of best fit for the 15Hz and 90Hz profiles for the VEXAT model (listed on Fig. 6A. and D.) and the Hill model (listed on Fig. 6G. and J., respectively).
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig6.tif/full/max/0/default.jpg)
The 15 Hz perturbations show that the VEXAT model’s performance is mixed: in the time-domain (A.) the VAF is lower than the 78-99% analyzed by Kirsch et al. [5]; the gain response (B.) follows the profile in Figure 3 of Kirsch et al. [5], while the phase response differs (C.). The response of the VEXAT model to the 90 Hz perturbations is much better: a VAF of 91% is reached in the time-domain (D.), the gain response follows the response of the cat soleus analyzed by Kirsch et al. [5], while the phase-response follows biological muscle closely for frequencies higher than 30 Hz. Although the Hill’s time-domain response to the 15 Hz signal has a higher VAF than the VEXAT model (G.), the RMSE of the Hill model’s gain response (H.) and phase response (I.) shows it to be more in error than the VEXAT model. While the VEXAT model’s response improved in response to the 90 Hz perturbation, the Hill model’s response does not: the VAF of the time-domain response remains low (J.), neither the gain (K.) nor phase responses (L.) follow the data of Kirsch et al. [5]. Note that the Hill model’s 90 Hz response was so nonlinear that the lowest frequency analyzed had to be raised from 4 Hz to 7 Hz to satisfy the criteria that (Cxy)2 ≥ 0.67.
The closeness of each model’s response to the spring-damper of best fit changes when a rigid-tendon is used instead of an elastic-tendon. While the VEXAT model’s response to the 15 Hz and 90 Hz perturbations improves slightly (compare Fig. 6A-F to Fig. 16A-F in Appendix F), the response of the Hill model to the 15 Hz perturbation changes dramatically with the time-domain VAF rising from 55% to 85% (compare Fig. 6G-L to Fig. 16G-L in Appendix F). Although the Hill model’s VAF in response to the 15 Hz perturbation improved, the frequency response contains mixed results: the rigid-tendon Hill model’s gain response is better (Fig. 16H in Appendix F), while the phase response is worse in comparison to the elastic-tendon Hill model. While the rigid-tendon Hill model produces a better time-domain response to the 15 Hz perturbation than the elastic-tendon Hill model, this improvement has been made with a larger phase shift between force and length than biological muscle [5].
The gain and phase profiles of both models deviate from the spring-damper of best fit due to the presence of nonlinearities, even for small perturbations. Some of the VEXAT model’s nonlinearities in this experiment come from the tendon model (compare Fig. 6A-F to Fig. 16A-F in Appendix F), since the response of the VEXAT model with a rigid-tendon stays closer to the spring-damper of best fit. The Hill model’s nonlinearities originate from the underlying expressions for stiffness and damping of the Hill model, which are particularly nonlinear with a rigid-tendon model (Fig. 16G-L in Appendix F) The stiffness of a Hill model’s CE
is heavily influenced by the partial derivative of
also suffers from high degrees of nonlinearity for small perturbations about v M = 0 since the slope of
and damping
of the VEXAT’s CE do not change so drastically because these terms do not depend on the slope of the force-length relation, or the force-velocity relation (see Appendix B.4 for derivation).
By repeating the stochastic perturbation experiments across a range of isometric forces, Kirsch et al. [5] were able to show that the stiffness and damping of a muscle varies linearly with the active tension it develops (see Figure 12 of [5]). We have repeated our simulations of Kirsch et al.’s [5] experiments at ten nominal forces (spaced evenly between 2.5N and 11.5 N) and compared how the VEXAT model and the Hill model’s stiffness and damping coefficients compare to Figure 12 of Kirsch et al. [5] (Fig. 7). The stiffness and damping profile of the VEXAT model deviates a little from Kirsch et al.’s data [5] because XE’s dynamics at 35 Hz are still influenced by the Hill model embedded in Eqn. 16 (see Appendix B.4). Despite this, the VEXAT model develops similar stiffness and damping profile with either a viscoelastic-tendon (Fig. 7A & B) or a rigid-tendon (Fig. 7C & D). In contrast, when the Hill model is coupled with an elastic-tendon both its stiffness and damping are larger than Kirsch et al.’s data [5] at the higher tensions (Fig. 7A and B). This pattern changes when simulating a Hill model with a rigid-tendon: the model’s stiffness is slightly negative (Fig. 7C), while the model’s final damping coefficient is nearly three times the value measured by Kirsch et al. (Fig. 7D). Though a negative stiffness may seem surprising, Eqn. 26 shows a negative stiffness is possible at the nominal CE length of these simulations: just past
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig7.tif/full/max/0/default.jpg)
When coupled with an elastic-tendon, the stiffness (A.) and damping (B.) coefficients of best fit of both the VEXAT model and a Hill model increase with the tension developed by the MTU. However, both the stiffness and damping of the elastic-tendon Hill model are larger than Kirsch et al.’s coefficients (from Figure 12 of [5]), particularly at higher tensions. When coupled with rigid-tendon the stiffness (C.) and damping (D.) coefficients of the VEXAT model remain similar, as the values for
When the VAF of the VEXAT and Hill model is evaluated across a range of nominal tensions (ten values from 2.5 to 11.5N), frequencies (15 Hz, 35 Hz, and 90 Hz), amplitudes (0.4mm, 0.8mm, and 1.6mm), and tendon types (rigid and elastic) two things are clear: first, that the VEXAT model’s 64-100% VAF is close to the 78-99% VAF reported by Kirsch et al. [5] while the Hill model’s 28-95% VAF differs (Fig. 8); and second, that there are systematic vari-ations in VAF, stiffness, and damping across the different perturbation magnitudes and frequencies (see Tables 5 and 5 in Appendix E). Both models produce worse VAF values when coupled with an elastic-tendon (Fig. 8A, B, and C), though the Hill model is affected most: the mean VAF of the elastic-tendon Hill model is 67% lower than the mean VAF of the rigid-tendon model for the 0.4 mm 15 Hz perturbations (Fig. 8A). While the VEXAT model’s lowest VAF occurs in response to the low frequency perturbations (Fig. 8A) with both rigid and elastic-tendons, the Hill model’s lowest VAF varies with both tendon type and frequency: the rigid-tendon Hill model has its lowest VAF in response to the 1.6 mm 90 Hz perturbations (Fig. 8C) while the elastic-tendon Hill model has its lowest VAF in response to the 0.4 mm 15 Hz perturbations (Fig. 8A). It is unclear if biological muscle displays systematic shifts in VAF since Kirsch et al. [5] did not report the VAF of each trial.
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig8.tif/full/max/0/default.jpg)
Kirsch et al. [5] noted that the VAF of the spring-damper model of best fit captured between 78-99% across all experiments. We have repeated the perturbation experiments to evaluate the VAF across a range of conditions: two different tendon models, three perturbation bandwidths (15 Hz, 35 Hz, and 90 Hz), three perturbation magnitudes (0.4 mm, 0.8 mm and 1.6 mm), and ten nominal force values (spaced evenly between 2.5N and 11.5N). Each bar in the plot shows the mean VAF across all ten nominal force values, with the whiskers extending to the minimum and maximum value that occurred in each set. The mean VAF of the VEXAT model changes by up to 36% depending on the condition, with the lowest mean VAF occurring in response to the 1.6 mm 15 Hz perturbation with an elastic-tendon (A.), and the highest mean VAF occurring in response to the 90 Hz perturbations with the rigid-tendon (C.). In contrast, the mean VAF of the Hill model varies by up to 67% depending on the condition, with the lowest VAF occurring in the 15 Hz 0.4 mm trial with the elastic-tendon (A.), and the highest value VAF occurring in the 15 Hz 0.4 mm trial with the rigid-tendon (A.).
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig9.tif/full/max/0/default.jpg)
Herzog and Leonard [7] actively lengthened (A.) cat soleus muscles on the descending limb of the force-length curve (where
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig10.tif/full/max/0/default.jpg)
In the VEXAT model we consider two different versions of the force-length relation for each titin segment (A): a linear extrapolation, and a WLC model extrapolation. Leonard et al. [8] observed that active myofibrils continue to develop increasing amounts of tension beyond actin-myosin overlap (B, grey lines with ±1 standard deviation shown). When this experiment is replicated using the VEXAT model (B., blue & magenta lines) and a Hill model (C. red lines), only the VEXAT model with the linear extrapolated titin model is able to replicate the experiment with the titin-actin bond slipping off of the actin filament at
3.2 Active lengthening on the descending limb
We now turn our attention to the active lengthening insitu experiments of Herzog and Leonard [7]. During these experiments cat soleus muscles were actively lengthened by modest amounts
When Herzog and Leonard’s [7] active ramp-lengthening (Fig. 9A) experiment is simulated, both models produce a force transient initially (Fig. 9B), but for different reasons. The VEXAT model’s transient is created when the lumped crossbridge spring (the
After the initial force transient, the response of the two models diverges (Fig. 9B): the VEXAT model continues to develop increased tension as it is lengthened, while the Hill model’s tension drops before recovering. The VEXAT model’s continued increase in force is due to the titin model: when activated, a section of titin’s PEVK region remains approximately fixed to the actin element (Fig. 1C). As a result, the ℓ 2 element (composed of part of PEVK segment and the distal Ig segment) continues to stretch and generates higher forces than it would if the muscle were being passively stretched. While both the elastic and rigid-tendon versions of the VEXAT model produce the same stereotypical ramp-lengthening response (Fig. 9B), the rigid-tendon model develops slightly more tension because the strain of the MTU is solely borne by the CE.
In contrast, the Hill model develops less force during lengthening when it enters a small region of negative stiffness (Fig. 9B and C) because the passive-force-length curve is too compliant to compensate for the negative slope of the active force-length curve. Similarly, the damping coefficient of the Hill model drops substantially during lengthening (Fig. 9D). Equation 27 and Figure 2B shows the reason that damping drops during lengthening: d f V/dv M, the slope of the line in Fig. 2B, is quite large when the muscle is iso-metric and becomes quite small as the rate of lengthening increases.
After the ramp stretch is completed (at time 3.4 seconds in Fig. 9B), the tension developed by the cat soleus recovers slowly, following a profile that looks strikingly like a first-order decay. The large damping coefficient acting between the titin-actin bond slows the force recovery of the VEXAT model. We have tuned the value of
Once activation is allowed to return to zero, Herzog and Leonard’s data shows that the cat soleus continues to develop a tension that is ΔfB above passive levels (Fig. 9B for t > 8.5s). The force ΔfB is known as passive force enhancement, and is suspected to be caused by titin bind-ing to actin [76]. Since we model titin-actin forces using an activation-dependent damper, when activation goes to zero our titin model becomes unbound from actin. As such, both our model and a Hill model remain ΔfB below the experimental data of Herzog and Leonard (Fig. 9B) after lengthening and activation have ceased.
3.3 Active lengthening beyond actin-myosin overlap
One of the great challenges that remains is to decompose how much tension is developed by titin (Fig. 1C) separately from myosin (Fig. 1B) in an active sarcomere. Leonard et al.’s [8] active-lengthening experiment provides some insight into this force distribution problem because they recorded active forces both within and far beyond actin-myosin overlap. Leonard et al.’s [8] data shows that active force continues to develop linearly during lengthening, beyond actin-myosin overlap, until mechanical failure. When activated and lengthened, the myofibrils failed at a length of
Since Leonard et al.’s experiment [8] was performed on skinned rabbit myofibrils and not on whole muscle, both the VEXAT and Hill models had to be adjusted prior to simulation (see Appendix H for parameter values). To simulate a rabbit myofibril we created a force-length curve [75] consistent with the filament lengths of rabbit skeletal muscle [60] (1.12µm actin, 1.63µm myosin, and 0.07µm z-line width) and fit the force-length relations of the two titin segments to be consistent with the structure measured by Prado et al. [58] of rabbit psoas titin consisting of a 70%-30% mix of a 3300kD and a 3400kD titin isoform (see Appendix B.3 for fitting details and Appendix H for parameter values). Since this is a simulation of a fibril, we used a rigid-tendon of zero length (equivalent to ignoring the tendon), and set the pennation angle to zero.
As mentioned in Sec. 2, because this experiment includes extreme lengths, we consider two different force-length relations for each segment of titin (Fig. 10A): a linear extrapolation, and an extension that follows the WLC model. While both versions of the titin model are identical up to
The Hill model was similarly modified, with the pennation angle set to zero and coupled with a rigid-tendon of zero length. Since the Hill model lacks an ECM element the passive-force-length curve was instead fitted to match the passive forces produced in Leonard et al.’s data [8]. No adjustments were made to the active elements of the Hill model.
When the slow active stretch (0.1µm/sarcomere/s) of Leonard et al.’s experiment is simulated [8] only the VEXAT model with the linear titin element can match the experimental data of Leonard et al. [8] (Fig. 10B). The Hill model cannot produce active force for lengths greater than
The WLC titin model is not able to reach the extreme lengths observed by Leonard et al. [8]. The distal segment of the WLC titin model approaches its contour length early in the simulation and ensures that the the titin-actin bond is dragged off the end of the actin filament at
This simulation has also uncovered a surprising fact: the myofibrils in Leonard et al.’s [8] experiments do not fail at
3.4 Force-length and force-velocity
Although the active portion of the Hill model is embedded in Eqn. 16, it is not clear if the VEXAT model can still replicate Hill’s force-velocity experiments [9] and Gordon et al.’s [10] force-length experiments. Here we simulate both of these experiments using the cat soleus model that we have used for the simulations described in Sec. 3.1 and compare the results to the force-length and force-velocity curves that are used in the Hill model and in Eqn. 16 of the VEXAT model.
Hill’s force-velocity experiment [9] is simulated by activating the model, and then by changing its length to follow a shortening ramp and a lengthening ramp. During short-ening experiments, the CE shortens from
The VEXAT model produces forces that differ slightly from the f V that is embedded in Eqn. 16 while the Hill model reproduces the curve (Fig. 11). The maximum short-ening velocity of the VEXAT model is slightly weaker than the embedded curve due to the series viscoelasticity of the XE element. Although the model can be made to converge to the f V curve more rapidly by decreasing τ S this has the undesirable consequence of degrading the low-frequency response of the model during Kirsch et al.’s experiments [5] (particularly Fig. 6C., and F.). These small differences can be effectively removed by scaling
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig11.tif/full/max/0/default.jpg)
When Hill’s [9] force-velocity experiment is simulated (A.), the VEXAT model produces a force-velocity profile (blue dots) that approaches zero more rapidly during shortening than the embedded profile f V(·) (red lines). By scaling
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig12.tif/full/max/0/default.jpg)
When Gordon et al.’s [10] passive and active force-length experiments are simulated, the VEXAT model (blue dots) and the Hill model (red lines) produce slightly different force-length curves (A.) and force responses in the time-domain (B.). The VEXAT model produces a right shifted active force-length curve, when compared to the Hill model due to the series elasticity of the XE element. By shifting the underlying curve by
Gordon et al.’s [10] force-length experiments were simulated by first passively lengthening the CE, and next by measuring the active force developed by the CE at a series of fixed lengths. Prior to activation, the passive CE was simulated for a brief period of time in a passive state to reduce any history effects due to the active titin element. To be consistent with Gordon et al.’s [10] experiment, we subtracted off the passive force from the active force before producing the active-force-length profile.
The simulation of Gordon et al.’s [10] experiment shows that the VEXAT model (Fig. 12A, blue dots) produces a force-length profile that is shifted to the right of the Hill model (Fig. 12A, red line) due to the series elasticity in-troduced by the XE. We can solve for the size of this right-wards shift by noting that Eqn. 16 will drive the
allowing us to solve for
the isometric strain of the XE. Since there are two viscoelastic XE elements per CE, the VEXAT model has an active force-length characteristic that shifted to the right of the embedded f L curve by a constant
4 Discussion
A muscle model is defined by the experiments it can replicate and the mechanisms it embodies. We have developed the VEXAT muscle model to replicate the force response of muscle to a wide variety of perturbations [5], [7], [8] while also retaining the ability to reproduce Hill’s force-velocity [9] experiment and Gordon et al.’s [10] force-length experiments. The model we have developed uses two mechanisms to capture the force response of muscle over a large variety of time and length scales: first, a viscoelastic crossbridge element that over brief time-scales appears as a springdamper, and at longer time-scales mimics a Hill-model; second, a titin element that is capable of developing active force during large stretches.
The viscoelastic crossbridge and titin elements we have developed introduce a number of assumptions into the model. While there is evidence that the activationdependent stiffness of muscle originates primarily from the stiffness of the attached crossbridges [62], the origins of the activation-dependent damping observed by Kirsch et al. [5] have not yet been established. We assumed that, since the damping observed by Kirsch et al. [5] varies linearly with activation, the damping originates from the attached crossbridges. Whether this damping is intrinsic or is due to some other factor remains to be established. Next, we have also assumed that the force developed by the XE converges to a Hill model [17] given enough time (Eqn. 16). A recent experiment of Tomalka et al. [80] suggests the force developed by the XE might decrease during lengthening rather than increasing as is typical of a Hill model [17]. If Tomalka et al.’s [80] observations can be replicated, the VEXAT model will need to be adjusted so that the the XE element develops less force during active lengthening while the active-titin element develops more force. Finally, we have assumed that actin-myosin sliding acceleration (due to crossbridge cycling) occurs when there is a force imbalance between the external force applied to the XE and the internal force developed by the XE as shown in Eqn. 16. This assumption is a departure from previous models: Hill-type models [16], [17] assume that the tension applied to the muscle instantaneously affects the actin-myosin sliding velocity; Huxley models [11] assume that the actin-myosin sliding velocity directly affects the rate of attachment and detachment of crossbridges.
The active titin model that we have developed makes assumptions similar to Rode et al. [38] and Schappacher-Tilp et al. [40]: some parts of the PEVK segment bond to actin, and this bond cannot do any positive work on titin. The assumption that the bond between titin and actin cannot do positive work means that titin cannot be actively preloaded: it can only develop force when it is stretched. In contrast, two mechanisms have been proposed that make it possible for titin to be preloaded by crossbridge cycling: Nishikawa’s [39] winding filament theory and DuVall et al.’s [44] titin entanglement hypothesis. If titin were significantly preloaded by crossbridge cycling, the titin load path would support higher forces and the myosin-actin load path would bear less force. Accordingly, the overall stiffness of the CE would be reduced, affecting our simulations of Kirsch et al. [5]: lower myosin-actin loads mean fewer attached crossbridges, since crossbridges are stiff in comparison to titin, the stiffness of the CE would decrease (see Appendix A). Hopefully experimental work will clarify if titin can be actively preloaded by crossbridges in the future.
Both the viscoelastic crossbridge and active titin elements include simple myosin-actin and titin-actin bond models that improve accuracy but have limitations. First, the viscoelastic crossbridge element has been made to represent a population of crossbridges in which the contribution of any single crossbridge is negligible. Though it may be possible for the XE model to accurately simulate a maximally activated single sarcomere (which has roughly 20 attached crossbridges per half sarcomere [12], [81]) the accuracy of the model will degrade as the number of attached cross-bridges decreases. When only a single crossbridge remains, the XE model’s output will be inaccurate because it can only generate force continuously while a real crossbridge generates force discretely each time it attaches to, and detaches from, actin. Next, we have used two equations, Eqns. 16 and 21, that assume myosin-actin and titin-actin interactions are temperature-invariant and scale linearly with size (
While the cat soleus XE and titin model parameters (Table 1 G, H, and I) can be used as rough default values, these parameters should be refit to accurately simulate muscle that differs in scale or temperature from cat soleus. Finally, the VEXAT model in its current form ignores phenomena related to submaximal contractions: the shift in the peak of the force-length relation [84], and the scaling of the maximum shortening velocity [85]. We hope to include these phenomena in a later version of the VEXAT model to more accurately simulate submaximal contractions.
The model we have proposed can replicate phenomena that occur at a wide variety of time and length scales: Kirsch et al.’s experiments [5] which occur over small time and length scales; and the active lengthening experiments of Herzog and Leonard [7] and Leonard et al. [8] which occur over physiological and supra-physiological length scales. In contrast, we have shown in Sec. 3.1 to 3.3 that a Hill-type model compares poorly to biological muscle when the same set of experiments are simulated. We expect that a Huxley model [11] is also likely to have difficulty reproducing Kirsch et al.’s experiment [5] because the model lacks an active damping element. Since titin was discovered [23] long after Huxley’s model was proposed [11], a Huxley model will be unable to replicate any experiment that is strongly influenced by titin such as Leonard et al.’s experiment [8].
Although there have been several more recent muscle model formulations proposed, none have the properties to simultaneously reproduce the experiments of Kirsch et al. [5], Herzog and Leonard [7], Leonard et al. [8], Hill [9], and Gordon et al. [10]. Linearized impedance models [14], [15] can reproduce Kirsch et al.’s experiments [5], but these models lack the nonlinear components needed to reproduce Gordon et al.’s force-length curve [10] and Hill’s forcevelocity curve [9]. The models of Forcinito et al. [18], and Tahir et al. [41] have a structure that places a contractile element in series with an elastic-tendon. While this is a commonly used structure, at high frequencies the lack of damping in the tendon will drive the phase shift between length and force to approach zero. The measurements and model of Kirsch et al. [5], in contrast, indicate that the phase shift between length and force approaches ninety degrees with increasing frequencies. Though the Hill-type models of Haeufle et al. [20] and Günther et al. [22] have viscoelastic tendons, these models have no representation, of the viscoelasticity of the CE’s attached crossbridges. Similar to the Hill-type muscle model evaluated in this work [17], it is likely that models of Haeufle et al. [20] and Günther et al. [22] will not be able to match the frequency response of biological muscle. De Groote et al. [51], [52] introduced a short-range-stiffness element in parallel to a Hill model to capture the stiffness of biological muscle. While De Groote et al.’s [51], [52] formulation improves upon a Hill model it is unlikely to reproduce Kirsch et al.’s experiment [5] because we have shown in Sec. 3.1 that a Hill model has a frequency response that differs from biological muscle. Rode et al.’s [38] muscle model also uses a Hill model for the CE and so we expect that this model will have the same difficulties reproducing Kirsch et al.’s [5] experiment. Schappacher-Tilp et al.’s model [40] extends a Huxley model [11] by adding a detailed titin element. Similar to a Huxley model, Schappacher-Tilp et al.’s model [40] will likely have difficulty reproducing Kirsch et al.’s experiment [5] because it is missing an active damping element.
While developing this model, we have come across open questions that we hope can be addressed in the future. How do muscle stiffness and damping change across the forcelength curve? Does stiffness and damping change with velocity? What are the physical origins of the active damping observed by Kirsch et al. [5]? What are the conditions that affect passive-force enhancement, and its release? In addition to pursuing these questions, we hope that other researchers continue to contribute experiments that are amenable to simulation, and to develop musculotendon models that overcome the limitations of our model. To help others build upon our work, we have made the source code of the model and all simulations presented in this paper available online14.
Acknowledgements
Financial support is gratefully acknowledged from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy (EXC 2075 – 390740016) through the Stuttgart Center for Simulation Science (SimTech), from DFG grant no. MI 2109/1-1, the Lighthouse Initiative Geriatronics by StMWi Bayern (Project X, grant no. 5140951), and NSERC of Canada.
A The stiffness of the actin-myosin and titin load paths
A single half-myosin can connect to the surrounding six actin filaments through 97.9 crossbridges. A 0.800 µm half-myosin has a pair crossbridges over 0.700 µm of its length every 14.3nm which amounts to 97.9 per halfmyosin [12]. Although 97.9 crossbridges does not make physical sense, here we will evaluate the stiffness of the CE assuming that fractional crossbridges exist and that attached crossbridges can be perfectly distributed among the 6 available actin filaments: the alternative calculation is more complicated and produces stiffness values that differ only in the 3rd significant digit. Assuming a duty cycle of 20% [81] (values between 5-90% have been reported [86]), at full actin-myosin overlap there will be 19.6 crossbridges attached to the 6 surrounding actin filaments. Assuming that these 19.6 crossbridges are evenly distributed between the 6 actin filaments, each single actin will be attached to 3.26 attached crossbridges.
At full overlap, the Z-line is 1 actin filament length LA (1.12 µm in rabbits [60]) from the M-line. The average point of crossbridge attachment is in the middle of the halfmyosin at a distance of 0.45 µm from the M-line (0.1 µm is bare and 0.35 µm is half of the remaining length), which is LA−0.45 µm from the Z-line. A single actin filament has a stiffness of 46 − 68 pN/nm [60] while a single crossbridge has a stiffness of 0.69 ± 0.47 pN/nm [62]. Since stiffness scales inversely with length, actin’s stiffness between the Zline and the average attachment point is 81.8−121 pN/nm. Finally, the stiffness of each actin filament and its 3.26 attached crossbridges is 0.712 − 3.67 pN/nm and all 6 together have a stiffness of 4.27 − 22.0 pN/nm.
Myosin has a similar stiffness as a single actin filament [61], with the section between the average attach-ment point and the M-line having a stiffness of 76.9 − 113 pN/nm. The final active stiffness of half-sarcomere is 4.05 − 18.4 pN/nm which comes from comes from the se-ries connection of the group of 6 actin filaments, with 19.6 crossbridges, and finally the single myosin filament. When this procedure is repeated assuming that only a single crossbridge is attached the stiffness drops to 0.22−1.15 pN/nm, which is slightly less than the stiffness of a single crossbridge 15.
The force-length profile of a single rabbit titin has been measured by Kellermayer et al. [87] using laser tweezers to apply cyclical stretches. By digitizing Fig. 4B (blue line) of Kellermayer et al. [87] we arrive at a stiffness for titin of 0.0058 − 0.0288 pN/nm at 2 µm (for a total sarcomere length of 4 µm or
When activated, the stiffness of our rabbit psoas linear-titin model (described in Sec. 3.3, fitted in Appendix B.3, and with the parameters shown in Appendix H) doubles, which would increase titin’s stiffness to 0.0696 − 0.346 pN/nm at 2 µm and 0.606 − 1.11 pN/nm at 4 µm.
Comparing the actin-myosin and titin stiffness ranges (Fig. 14) makes it clear that the stiffness of actin-myosin with 1 attached crossbridge (AM:Low in Fig. 14) is comparable to the highest stiffness values we have estimated for titin (TA:High in Fig. 14). When all 20% of the available crossbridges are attached (AM:High in Fig. 14), the average stiffness of the actin-myosin load path is roughly one order of magnitude stiffer than the highest stiffness values of titin (TA:High in Fig. 14), and two to three orders of magnitude higher than the lowest stiffness titin load path (TP:Low in Fig. 14). Similarly, the maximum XE stiffness and titin stiffness in this work are separated by roughly an order of magnitude: the cat soleus model has a XE stiff-ness of
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig13.tif/full/max/0/default.jpg)
To evaluate the stiffness of the actin-myosin load path, we first determine the average point of attachment. Since the actin filament length varies across species we label it LA. Across rabbits, cats and human skeletal muscle myosin geometry is consistent [75]: a half-myosin is 0.8µm in length with a 0.1µm bare patch in the middle. Thus at full overlap the average point of attachment is 0.45µm from the M-line, or LA−0.45µm from the Z-line at
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig14.tif/full/max/0/default.jpg)
The stiffness of a rabbit’s actin-myosin load path with a single attached crossbridge (1 XB) exceeds the stiff-ness of its titin filament at lengths of
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig15.tif/full/max/0/default.jpg)
The normalized tendon force-length curve (A) has been been fit to match the cat soleus tendon stiffness measurements of Scott and Loeb [74]. The data of Netti et al. [72] allow us to develop a model of tendon damping as a linear function of tendon stiffness. By normalizing the measurements of Netti et al. [72] by the maximum storage modulus we obtain curves that are equivalent to the normalized stiffness (B) and damping (C) of an Achilles tendon from a rabbit. Both normalized tendon stiffness and damping follow similar curves, but at different scales, allowing us to model tendon damping as a linear function of tendon stiffness (C).
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig16.tif/full/max/0/default.jpg)
When coupled with a rigid-tendon, the VEXAT model’s VAF (A.), gain response (B.), and phase response (C.) more closely follows the data of Kirsch et al. (Figure 3) [5] than when an elastic-tendon is used. This improvement in accuracy is also observed at the 90 Hz perturbation (D., E., and F.), though the phase response of the model departs from Kirsch et al.’s data [5] for frequencies lower than 30 Hz. Parts of the Hill model’s response to the 15 Hz perturbation are better with a rigid-tendon, with a higher VAF (G.), a lower RMSE gain-response (H.). but have a poor phase-response (I.). In response to the higher frequency perturbations, the Hill model’s response is poor with an elastic (see Fig. 6) or rigid-tendon. The VAF in response to the 90 Hz perturbation remains low (J.), and neither the gain (K.) nor the phase response of the Hill model (L.) follow the data of Kirsch et al. [5]. The rigid-tendon Hill model’s nonlinearity was so strong that the lowest frequency analyzed had to be raised from 4 Hz to 21 Hz to meet the criteria that (Cxy)2 ≥ 0.67.
B Model Fitting
Many of the experiments simulated in this work [5], [7] have been performed using cat soleus muscle. While we have been able to extract some architectural parameters directly from the experiments we simulate (
B.1 Fitting the tendon’s stiffness and damping
Similar to previous work [17], we model the force-length relation of the tendon using a quintic Bézier spline (Fig. 15A) that begins at
The similarity in shape is likely not a coincidence.
The nonlinear characteristics (Fig. 15) tendon originates from its microscopic structure. Tendon is composed of many fiber bundles with differing resting lengths [72]. Initially the tendon’s fiber bundles begin crimped, but gradually stretch as the tendon lengthens, until finally all fiber bundles are stretched and the tendon achieves its maximum stiffness (Fig. 15B) and damping (Fig. 15C) [72]. Accordingly, in Eqn. 23 we have described the normalized damping of the tendon as being equal to the normalized stiffness of the tendon scaled by a constant U. To estimate U we have used the measurements of Netti et al. [72] (Fig. 15 B and C) and solved a least-squares problem
to arrive at a value of U = 0.057. The resulting damping model (Fig. 15C) fits the measurements of Netti et al. [72] closely.
B.2 Fitting the CE’s Impedance
We can now calculate the normalized impedance of the XE using the viscoelastic-tendon model we have constructed and Kirsch et al.’s [5] measurements of the impedance of the entire MTU. Since an MTU is structured with a CE in series with a tendon, the compliance of the MTU is given by
where
we can use a similar procedure to evaluate
Using Eqn. 34, we arrive at
The stiffness of the CE along the tendon is
which can be expanded to
Since we are using a constant thickness pennation model
and thus
which simplifies to
Similarly, the constant thickness pennation model means that
which leads to
Recognizing that
we can solve for k M in terms of
We can use a similar process to transform
which expands to a much smaller expression
than Eqn. 36 since α does not depend on v M, and thus ∂α/∂v M = 0. By taking a time derivative of Eqn. 40 we arrive at
which allows us to solve for
By recognizing that
and using Eqns. 44 and 46 we can evaluate β M in terms of
B.3 Fitting the force-length curves of titin’s segments
The nonlinear force-length curves used to describe titin (
First, we fit the passive force-length curve to the data of Herzog and Leonard [7] to serve as a reference. The curve f PE begins at the normalized length and force coordinates of
the squared differences between f PE and the passive forcelength data of Herzog and Leonard [7] (Fig. 2A shows both the data and the fitted f PE curve) are minimized. While f PE is not used directly in the model, it serves as a useful reference for constructing the ECM and titin force-length curves. We assume that the ECM force-length curve is a linear scaling of f PE
where P is a constant. In this work, we set P to 56% which is the average ECM contribution that Prado et al. [58] measured across 5 different rabbit skeletal muscles16. The remaining fraction, 1 − P, of the force-length curve comes from titin.
In mammalian skeletal muscle, titin has three elastic segments [58] connected in series: the proximal Ig segment, the PEVK segment, and the distal Ig segment that is between the PEVK region and the myosin filament (Fig. 1A). Trombitás et al. [26] labelled the PEVK region of titin with antibodies allowing them to measure the distance between the Z-line and the proximal Ig/PEVK boundary (Zℓ IgP/PEVK), and the distance between the Z-line and the PEVK/distal Ig boundary (Zℓ PEVK/IgD), while the passive sarcomere was stretched from 2.35 − 4.46µm. By fitting functions to Trombitás et al.’s [26] data we can predict the length of any of titin’s segments under the following assumptions: the T12 segment is rigid (Fig. 1A), the distal Ig segment that overlaps with myosin is rigid (Fig. 1A), and that during passive stretching the tension throughout the titin filament is uniform. Since the sarcomeres in Trombitás et al.’s [26] experiments were passively stretched it is reasonable to assume that tension throughout the free part of the titin filament is uniform because the bond between titin and actin depends on calcium [31], [36] and crossbridge attachment [8].
We begin by digitizing the data of Figure 5 of Trombitás et al. [26] and using the least-squares method to fit lines to Zℓ IgP/PEVK and Zℓ PEVK/IgD (where the superscripts mean fromℓto and so Zℓ IgP/PEVK is the distance from the Z-line to the border of the IgP/PEVK segments). From these lines of best fit we can evaluate the normalized length of the proximal Ig segment
the normalized length of the PEVK segment
and the normalized length of the distal Ig segment
as a function of sarcomere length. Next, we extract the coefficients for linear functions that evaluate the lengths of
given the
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl2.tif/full/max/0/default.jpg)
The coefficients of the normalized lengths of
These functions can be scaled to fit a titin filament of a differing geometry. Many of the experiments simulated in this work used cat soleus. Although the lengths of titin filament segments in cat soleus have not been measured, we assume that it is a scaled version of a human soleus titin filament (68 proximal Ig domains, 2174 PEVK residues, and 22 distal Ig domains [26]) since both muscles contain predominately slow-twitch fibers: slow twitch fibers tend to express longer, more compliant titin filaments [58]. Since the optimal sarcomere length in cat skeletal muscle is shorter than in human skeletal muscle (2.43 µm vs. 2.73 µm, [75]) the coefficients for Eqns. 53-55 differ slightly (see the feline soleus column in Table 2). In addition, by assuming that the titin filament of cat skeletal muscle is a scaled version of the titin filament found in human skeletal muscle, we have implicitly assumed that the cat’s skeletal muscle titin filament has 60.5 proximal Ig domains, 1934.7 PEVK residues, and 19.6 distal Ig domains. Although a fraction of a domain does not make physical sense, we have not rounded to the nearest domain and residue to preserve the sarcomere length-based scaling.
In contrast, the rabbit psoas fibril used in the simulation of Leonard et al. [8] has a known titin geometry (50 proximal Ig domains, 800 PEVK residues, and 22 distal Ig domains [58]) which differs substantially from the isoform of titin expressed in the human soleus. To create the rabbit psoas titin length functions
is what remains from the half-sarcomere once the rigid lengths of titin (0.100 µm for LT12 and 0.8150 µm for LIgD [60]) and the PEVK segment length have been subtracted away. The function that describes
which produce the coefficients that appear in the rabbit psoas column in Table 2. While we have applied this approach to extend Trombitás et al.’s [26] results to a rabbit psoas, in principle this approach can be applied to any isoform of titin provided that its geometry is known, and the Ig domains and PEVK residues in the target titin behave similarly to those in human soleus titin.
The only detail that remains is to establish the shape of the IgP, PEVK, and IgD force-length curves. Studies of individual titin filaments, and of its segments, make it clear that titin is mechanically complex. For example, the tandem Ig segments (the IgD and IgP segments) are composed of many folded domains (titin from human soleus has two Ig segments that together have nearly 100 domains [26]). Each domain appears to be a simple nonlinear spring until it unfolds and elongates by nearly 25 nm in the process [88]. Unfolding events appear to happen individually during lengthening experiments [88], with each unfolding event occurring at a slightly greater tension than the last, giving an Ig segment a force-length curve that is saw-toothed. Although detailed models of titin exist that can simulate the folding and unfolding of individual Ig domains, this level of detail comes at a cost of a state for each Ig domain which can add up to nearly a hundred extra states [40] in total.
Active and passive lengthening experiments at the sarcomere-level hide the complexity that is apparent when studying individual titin filaments. The experiments of Leonard et al. [8] show that the sarcomeres in a filament (from rabbit psoas) mechanically fail when stretched passively to an average length of
Since we are interested in a computationally efficient model that is accurate at the whole muscle level, we model titin as a multi-segmented nonlinear spring but omit the states needed to simulate the folding and unfolding of Ig domains. Simulations of active lengthening using our titin model will exhibit the enhanced force development that appears in experiments [7], [8], but will lack the nonlinear saw-tooth force-length profile that is measured when individual titin filaments are lengthened [88]. To have the broadest possible application, we will fit titin’s force-length curves to provide reasonable results for both moderate [7] and large active stretches [8]. Depending on the application, it may be more appropriate to use a stiffer force-length curve for the Ig segment if normalized sarcomere lengths stays within
To ensure that the serially connected force-length curves of
which is formed by the series connection of
Since each of titin’s segments is exposed to the same tension in Trombitás et al.’s experiment [26] the slopes of the lines that Eqns. 53-55 describe are directly proportional to the relative compliance (inverse of stiffness) between of each of titin’s segments. Using this fact, we can define the normalized stretch rates of the proximal titin segment
and the distal titin segment
which are proportional to the compliance of two titin segments in our model. When both the
and
Dividing Eqn. 63 by 64 eliminates the unknown
to the relative changes in Eqns. 61 and 62. Substituting Eqns. 65, and 60 into Eqn. 59 yields the expression
which can be simplified to
and this expression can be evaluated using the terminal stiffness of titin
The curves
and ℓ 2 segment
at a CE length of
By construction, the spring network formed by the
The process we have used to fit the ECM and titin’s segments makes use of data within modest normalized CE lengths (2.35-4.46µm, or
We have modified the WLC to include a slack length
where B is a scaling factor and the normalized segment length
where
by counting the number of proximal Ig domains (N IgP), at the number of PEVK residues (QN PEVK) associated with ℓ 1 and by scaling each by the maximum contour length of each Ig domain (25nm [88]), and each PEVK residue (between 0.32 [65] and 0.38 nm [90] see pg. 254). This contour length defines the maximum length of the segment, when all of the Ig domains and PEVK residues have been unfolded. Similarly, the contour length of
Next, we define the slack length by linearly extrapolating backwards from the final fitted force (1 − P)
and similarly
We can now solve for B in Eqn. 71 so that
B.4 Fitting the XE’s Impedance
With the passive curves established, we can return to the problem of identifying the normalized maximum stiffness
By noting that all of our chosen state variables in Eqn. 13 are independent and by making use of the kinematic relationships in Eqns. 9 and 10 we can reduce Eqn. 77 to
and solve for
When using to the data from Figure 12 in Kirsch et al. [5], we end up with
As with k M, the expression for
which evaluates to
The dimensionless parameters
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl3.tif/full/max/0/default.jpg)
Normalized titin and crossbridge parameters fit to data from the literature.
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl4.tif/full/max/0/default.jpg)
Mean normalized stiffness coefficients (A.), mean normalized damping coefficients (B.), VAF (C.), and the bandwidth (D.) of linearity (coherence squared > 0.67) for models with elastic-tendons. Here the proposed model has been fitted to Figure 12 of Kirsch et al. [5], while the experimental data from Kirsch et al. [5] comes from Figures 9 and 10. Experimental data from Figure 12 from Kirsch et al. has not been included in this table because it would only contribute 1 entry and would overwrite values from Figures 9 and 10. The impedance experiments at each combination of perturbation amplitude and frequency have been evaluated at 10 different nominal forces linearly spaced between 2.5N and 11.5N. The results presented in the table are the mean values of these ten simulations. The VAF is evaluated between the model and the spring-damper of best fit, rather than to the response of biological muscle (which was not published by Kirsch et al. [5]). Finally, model values for the VAF (C.) and the bandwidth of linearity (D.) that are worse than those published by Kirsch et al. [5] appear in bold font.
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl5.tif/full/max/0/default.jpg)
Mean normalized stiffness coefficients (A.), mean normalized damping coefficients (B.), VAF (C.), and the bandwidth (D.) of linearity (coherence squared > 0.67) for models with rigid tendons. All additional details are identical to those of Table except the tendon of the model is rigid.
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_tbl6.tif/full/max/0/default.jpg)
The VEXAT and Hill model’s fitted rabbit psoas fibril MTU parameters. As in Table 6, parameters shared by the VEXAT and Hill model are highlighted in grey. Short forms are used to save space: length ‘len’, velocity ‘vel’, acceleration ‘acc’, half ‘h’, activation ‘act’, segment ‘seg’, threshold ‘thr’, and stiffness ‘stiff’. The letter preceding a reference indicates the experimental animal:’C’ for cat, ‘H’ for human, while nothing at all is rabbit skeletal muscle. Letters following a reference indicate how the data was used to evaluate the parameter: ‘A’ for arbitrary for simulating Leonard et al. [8], ‘n/a’ for a parameter that is not applicable to a fibril model, ‘—’ value taken from the cat soleus MTU, ‘C’ calculated, ‘F’ fit, ‘E’ estimated, ‘S’ scaled, and ‘D’ for default if a default value from another model was used. Only parameters that do not affect the outcome of our simulation of Leonard et al. [8] are marked ‘A’. Clearly the parameters that appear in this Table do not represent a generic rabbit psoas fibril model, but instead a rabbit psoas fibril model that is sufficient to simulate the experiment of Leonard et al. [8]. Finally, values for N IgP, N PEVK, and N IgD were obtained by taking a 70% and 30% average of the values for 3300 kD and 3400 kD titin to match the composition of rabbit psoas titin as closely as possible.
C Model Initialization
Solving for an initial state is challenging since we are given a, ℓ P, and v P and must solve for v S, ℓ S, and ℓ 1 for a rigid-tendon model, and additionally ℓ M if an elastic-tendon model is used. The type of solution that we look for is one that produces no force or state transients soon after a simulation begins in which activation and path velocity is well approximated as constant. Our preliminary simulations found that satisfactory solutions were found by iterating over both
In the outer loop, we iterate over values of
In the inner loop, we iterate over values of
and we use fixed-point iteration to solve for
D Evaluating a muscle model’s frequency response
To analyze the the frequency response of a muscle to length perturbation we begin by evaluating the length change
and force change
with respect to the nominal length
where ∗ is the convolution operator. Each of these signals can be transformed into the frequency-domain [78] by taking the Fourier transform ℱ (·) of the time-domain signal, which produces a complex (with real and imaginary parts) signal. Since convolution in the time-domain corresponds to multiplication in the frequency-domain, we have
In Eqn. 87 we are interested in solving for H(s). While it might be tempting to evaluate H(s) as
the result will poorly estimate H(s) because Y (s) is only approximated by H(s) X(s): Y (s) may contain nonlinearities, non-stationary signals, and noise that cannot be described by H(s) X(s).
Using cross-spectral densities, Koopmans [79] (p. 140) derived the estimator
that minimizes the squared errors between Y (s) and its linear approximation of H(s) X(s). The cross-spectral density Gxy between x(t) and y(t) is given by
the Fourier transform of the cross-correlation (⋆) between x(t) and y(t). When the order of x(t) and y(t) are reversed in Eqn. 90 the result is Gyx, while Gxx and Gyy are produced by taking the Fourier transform of x(t) ⋆ x(t) and y(t) ⋆ y(t) respectively.
Though Koopmans’s [79] estimator is a great improvement over Eqn. 88, the accuracy of the estimate can be further improved using Welch’s method [91]. Welch’s method [91] breaks up the time domain signal into K segments, transforms each segment into the frequency domain, and returns the average across all segments. Using Welch’s method [91] with K segments allows us to evaluate
which has a lower frequency resolution than Eqn. 89, but an improved accuracy in H(s). Now we can evaluate the gain of H(s) as
while the phase of H(s) is given by
where ℝ (H(s)) and 𝕀 (H(s)) are the real and imaginary parts of H(s) respectively.
The transfer function estimated in Eqn. 91 is meaningful only when y(t) can be approximated as a linear time-invariant function of x(t). By evaluating the coherence [79] (p. 137) between x(t) and y(t)
we can determine the strength of the linear association between X(s) and Y (s) at each frequency. When Cxy is close to 1 it means that Y (s) is well approximated by H(s) X(s). As Cxy approaches 0, it means that the approximation of Y (s) by H(s) X(s) becomes poor.
Kirsch et al. [5] analyzed a bandwidth that spanned from 4 Hz up to the cutoff frequency of the low-pass filter ap-plied to the input signal x(t) (15 Hz, 35 Hz, and 90 Hz).
Unfortunately, we cannot use this bandwidth directly when analyzing model output because we have no guarantee that the simulated output is sufficiently linear in this range. Instead, to strike a balance between accuracy and consistency with Kirsch et al. [5], we analyze the bandwidth that is common to Kirsch et al.’s [5] defined range and has the minimum acceptable (Cxy)2 of 0.67 that is pictured in Fig. 3 of Kirsch et al.
E Simulation summary data of Kirsch et al
F Supplementary plots: Gain and phase response rigid-tendon muscle models
G Supplementary plots: active lengthening on the descending limb
H Rabbit psoas fibril model parameters
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig17.tif/full/max/0/default.jpg)
Simulation results of the 3 mm/s (A.) active lengthening experiment of Herzog and Leonard [7] (B.). As with the 9 mm/s trial, the Hill model’s force response drops during the ramp due to a small region of negative stiffness introduced by the descending limb of the force-length curve (C.), and a reduction in damping (D.) due to the flattening of the force-velocity curve. Note: neither model was fitted to this trial.
![](https://prod--epp.elifesciences.org/iiif/2/88344%2Fv3%2Fcontent%2F534347v6_fig18.tif/full/max/0/default.jpg)
Simulation results of the 27 mm/s (A.) active lengthening experiment of Herzog and Leonard [7] (B.). As with the prior simulations the Hill model exhibits a small region of negative stiffness introduced by the descending limb of the force-length curve (C.) and a drop in damping (D.). Note: neither model was fitted to this trial.
References
- [1]Adaptation to stable and unstable dynamics achieved by combined impedance control and inverse dynamics modelJournal of neurophysiology 90:3270–3282https://doi.org/10.1152/jn.01112.2002
- [2]Use of self-selected postures to regulate multijoint stiffness during unconstrained tasksPloS one 4https://doi.org/10.1371/journal.pone.0005411
- [3]Impedance control reduces instability that arises from motor noiseJournal of neuroscience 29:12–12https://doi.org/10.1523/JNEUROSCI.2826-09.2009
- [4]The central nervous system stabilizes unstable dynamics by learning optimal impedanceNature 414:446–449https://doi.org/10.1038/35106566
- [5]Muscle stiffness during transient and continuous movements of cat muscle: Perturbation characteristics and physiological relevanceIEEE Transactions on Biomedical Engineering 41:758–770https://doi.org/10.1109/10.310091
- [6]The short range stiffness of active mammalian muscle and its effect on mechanical propertiesThe Journal of physiology 240:331–350https://doi.org/10.1113/jphysiol.1974.sp010613
- [7]Force enhancement following stretching of skeletal muscle: A new mechanismThe Journal of Experimental Biology 205:1275–1283https://doi.org/10.1242/jeb.205.9.1275
- [8]An activatable molecular spring reduces muscle tearing during extreme stretchingJournal of biomechanics 43:3063–3066https://doi.org/10.1016/j.jbiomech.2010.07.016
- [9]The heat of shortening and the dynamics constants of musclein Proceedings of the Royal Society of London 126:136–195https://doi.org/10.1098/rspb.1938.0050
- [10]The variation in isometric tension with sarcomere length in vertebrate muscle fibresThe Journal of Physiology 184:170–192https://doi.org/10.1113/jphysiol.1966.sp007909
- [11]Muscle structure and theories of contractionProgress in biophysics and biophysical chemistry 7:255–318https://doi.org/10.1016/S0096-4174(18)30128-8
- [12]Proposed mechanism of force generation in striated muscleNature 223:533–538https://doi.org/10.1038/233533a0
- [13]Fibersim: A flexible open-source model of myofilament-level contractionBiophysical Journal 121:175–182https://doi.org/10.1016/j.bpj.2021.12.021
- [14]The mechanics of multi-joint posture and movement controlBiological Cybernetics 52:315–331https://doi.org/10.1007/BF00355754
- [15]Neural, mechanical, and geometric factors subserving arm posture in humansJournal of Neuroscience 5:2732–2743https://doi.org/10.1523/JNEUROSCI.05-10-02732.1985
- [16]Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor controlCritical Reviews in Biomedical Engineering 17:359–411
- [17]Flexing computational muscle: Modeling and simulation of musculotendon dynamicsJournal of biomechanical engineering 135https://doi.org/10.1115/1.4023390
- [18]Can a rheological muscle model predict force depression/enhancement?Journal of biomechanics 31:1093–1099https://doi.org/10.1016/S0021-9290(98)00132-8
- [19]Architecture of the hind limb muscles of cats: Functional significanceJournal of Morphology 173:185–195https://doi.org/10.1002/jmor.1051730206
- [20]Hill-type muscle model with serial damping and eccentric force–velocity relationJournal of biomechanics 47:1531–1536https://doi.org/10.1016/j.jbiomech.2014.02.009
- [21]Is equilibrium point control feasible for fast goal-directed single-joint movements?Journal of Neurophysiology 95:2898–2912https://doi.org/10.1152/jn.00983.2005
- [22]The basic mechanical structure of the skeletal muscle machinery: One model for linking microscopic and macroscopic scalesJournal of theoretical biology 456:137–167https://doi.org/10.1016/j.jtbi.2018.07.023
- [23]Connectin, an elastic protein from myofibrilsThe Journal of Biochemistry 80:405–407
- [24]Titin: Major myofibrillar components of striated muscleProceedings of the National Academy of Sciences 76:3698–3702https://doi.org/10.1073/pnas.76.8.3698
- [25]Connectin filaments link thick filaments and z lines in frog skeletal muscle as revealed by immunoelectron microscopyThe Journal of cell biology 101:2167–2172https://doi.org/10.1083/jcb.101.6.2167
- [26]Pevk extension of human soleus muscle titin revealed by immunolabeling with the anti-titin antibody 9d10Journal of structural biology 122:188–196https://doi.org/10.1006/jsbi.1998.3984
- [27]Studies of the interaction between titin and myosinJournal of Cell Biology 131:1471–1481https://doi.org/10.1083/jcb.131.6.1471
- [28]Extensive eccentric contractions in intact cardiac trabeculae: Revealing compelling differences in contractile behaviour compared to skeletal musclesProceedings of the Royal Society B 286:1–9https://doi.org/10.1098/rspb.2019.0719
- [29]Residual and passive force enhancement in skinned cardiac fibre bundlesJournal of Biomechanics 109https://doi.org/10.1016/j.jbiomech.2020.109953
- [30]Huxleys’ missing filament: Form and function of titin in vertebrate striated muscleAnnual review of physiology 79:145–166https://doi.org/10.1146/annurev-physiol-022516-034152
- [31]Calcium-dependent inhibition of in vitro thin-filament motility by native titinFEBS letters 380:281–286https://doi.org/10.1016/0014-5793(96)00055-5
- [32]Binding of a native titin fragment to actin is regulated by pip2FEBS letters 429:95–98https://doi.org/10.1016/S0014-5793(98)00572-9
- [33]Regulation of the actin–myosin interaction by titinEuropean journal of biochemistry 271:4572–4581https://doi.org/10.1111/j.1432-1033.2004.04429.x
- [34]Interaction forces between f-actin and titin pevk domain measured with optical tweezersBiophysical Journal 93:2102–2109https://doi.org/10.1529/biophysj.107.106153
- [35]Differential actin binding along the pevk domain of skeletal muscle titinJournal of cell science 117:5781–5789https://doi.org/10.1242/jcs.01501
- [36]Calcium increases titin n2a binding to f-actin and regulated thin filamentsScientific reports 8:1–11https://doi.org/10.1038/s41598-018-32952-8
- [37]Differences in titin segmental elongation between passive and active stretch in skeletal muscleJournal of Experimental Biology 220:4418–4425https://doi.org/10.1242/jeb.160762
- [38]Titin-induced force enhancement and force depression: A ‘sticky-spring’mechanism in muscle contractions?Journal of theoretical biology 259:350–360https://doi.org/10.1016/j.jtbi.2009.03.015
- [39]Is titin a ‘winding filament’? a new twist on muscle contractionProceedings of the royal society B: Biological sciences 279:981–990https://doi.org/10.1098/rspb.2011.1304
- [40]A novel three-filament model of force generation in eccentric contraction of skeletal musclesPloS one 10https://doi.org/10.1371/journal.pone.0117634
- [41]Case study: A bio-inspired control algorithm for a robotic foot-ankle prosthesis provides adaptive control of level walking and stair ascentFrontiers in Robotics and AI 5https://doi.org/10.3389/frobt.2018.00036
- [42]A multi-scale continuum model of skeletal muscle mechanics predicting force enhancement based on actin–titin interactionBiomechanics and modeling in mechanobiology 15:1423–1437https://doi.org/10.1007/s10237-016-0772-7
- [43]A continuum-mechanical skeletal muscle model including actin-titin interaction predicts stable contractions on the descending limb of the force-length relationPLoS computational biology 13https://doi.org/10.1371/journal.pcbi.1005773
- [44]I-band titin interaction with myosin in the muscle sarcomere during eccentric contraction: The titin entanglement hypothesisBiophysical Journal 110https://doi.org/10.1016/j.bpj.2015.11.1625
- [45]Identification of intrinsic and reflex contributions to human ankle stiffness dynamicsIEEE transactions on biomedical engineering 44:493–504https://doi.org/10.1109/10.581944
- [46]Multijoint dynamics and postural stability of the human armExperimental brain research 157:507–517https://doi.org/10.1007/s00221-004-1864-7
- [47]NM-Clab, a model to assess the contributions of muscle visco-elasticity and afferent feedback to joint dynamicsJournal of biomechanics 41:1659–1667https://doi.org/10.1016/j.jbiomech.2008.03.014
- [48]A computational model of limb impedance control based on principles of internal model uncertaintyPloS one 5https://doi.org/10.1371/journal.pone.0013601
- [49]Modeling and simulating the neuromuscular mechanisms regulating ankle and knee joint stiffness during human locomotionJournal of neurophysiology 114:2509–2527https://doi.org/10.1152/jn.00989.2014
- [50]Bridging scales: A three-dimensional electromechanical finite element model of skeletal muscleSIAM Journal on Scientific Computing 30:2882–2904https://doi.org/10.1137/070691504
- [51]Contribution of muscle short-range stiffness to initial changes in joint kinetics and kinematics during perturbations to standing balance: A simulation studyJournal of biomechanics 55:71–77https://doi.org/10.1016/j.jbiomech.2017.02.008
- [52]Interaction between muscle tone, short-range stiffness and increased sensory feedback gains explains key kinematic features of the pendulum test in spastic cerebral palsy: A simulation studyPloS one 13https://doi.org/10.1371/journal.pone.0205763
- [53]Huxley-type cross-bridge models in largeish-scale musculoskeletal models; an evaluation of computational costJournal of biomechanics 83:43–48https://doi.org/10.1016/j.jbiomech.2018.11.021
- [54]In vivo specific tension of human skeletal muscleJournal of applied physiology 90:865–872https://doi.org/10.1152/jappl.2001.90.3.865
- [55]Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb musclesJournal of Biomechanics 44:109–115https://doi.org/10.1016/j.jbiomech.2010.08.033
- [56]Maximum velocity of shortening of three fibre types from horse soleus muscle: Implications for scaling with body sizeThe Journal of Physiology 431:173–185https://doi.org/10.1113/jphysiol.1990.sp018325
- [57]Are titin properties reflected in single myofibrils?Journal of biomechanics 45:1893–1899https://doi.org/10.1016/j.jbiomech.2012.05.021
- [58]Isoform diversity of giant proteins in relation to passive and active contractile properties of rabbit skeletal musclesThe Journal of general physiology 126:461–480https://doi.org/10.1085/jgp.200509364
- [59]Stiffness, working stroke, and force of single-myosin molecules in skeletal muscle: Elucidation of these mechanical properties via nonlinear elasticity evaluationCellular and molecular life sciences 70:4275–4292https://doi.org/10.1007/s00018-013-1353-x
- [60]Compliance of thin filaments in skinned fibers of rabbit skeletal muscleBiophysical journal 69:1000–1010https://doi.org/10.1016/S0006-3495(95)79975-1
- [61]X-ray evidence for the elongation of thin and thick filaments during isometric contraction of a molluscan smooth muscleJournal of Muscle Research & Cell Motility 15:659–671https://doi.org/10.1007/BF00121073
- [62]The stiffness of rabbit skeletal actomyosin cross-bridges determined with an optical tweezers transducerBiophysical journal 75:1424–1438https://doi.org/10.1016/S0006-3495(98)74061-5
- [63]Depression of cat soleus forces following isokinetic shorteningJournal of biomechanics 30:865–872https://doi.org/10.1016/S0021-9290(97)00046-8
- [64]Evaluation of direct collocation optimal control problem formulations for solving the muscle redundancy problemAnnals of biomedical engineering 44:2922–2936https://doi.org/10.1007/s10439-016-1591-9
- [65]Titin extensibility in situ: Entropic elasticity of permanently folded and permanently unfolded molecular segmentsThe Journal of cell biology 140:853–859https://doi.org/10.1083/jcb.140.4.853
- [66]Force enhancement in single skeletal muscle fibers on the ascending limb of the force-length relationshipThe Journal of Experimental Biology 207:2787–2791https://doi.org/10.1242/jeb.01095
- [67]Does residual force enhancement increase with increasing stretch magnitudes?Journal of biomechanics 42:1488–1492https://doi.org/10.1016/j.jbiomech.2009.03.046
- [68]Residual force enhancement is preserved for conditions of reduced contractile forceMedicine and science in sports and exercise 50:1186–1191https://doi.org/10.1249/mss.0000000000001563
- [69]The effect of phosphate and calcium on force generation in glycerinated rabbit skeletal muscle fibers. a steady-state and transient kinetic studyJournal of Biological Chemistry 265:20–20https://doi.org/10.1016/S0021-9258(17)30494-5
- [70]Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adultsJournal of Biomechanical Engineering 125:70–77https://doi.org/10.1115/1.1531112
- [71]Evidence that maximum muscle stress is not a constant: Differences in specific tension in elbow flexors and extensorsMedical engineering & physics 17:529–536https://doi.org/10.1016/1350-4533(95)00005-8
- [72]Structure-mechanical properties relationship of natural tendons and ligamentsJournal of Materials Science: Materials in Medicine 7:525–530https://doi.org/10.1007/BF00122175
- [73]Mechanics of feline soleus: I. effect of fascicle length and velocity on force outputJournal of Muscle Research & Cell Motility 17:207–219https://doi.org/10.1002/jmor.1052240109
- [74]Mechanical properties of aponeurosis and tendon of the cat soleus muscle during whole-muscle isometric contractionsJournal of Morphology 224:73–86https://doi.org/10.1002/jmor.1052240109
- [75]Length dependence of active force production in skeletal muscleJournal of applied physiology 86:1445–1457https://doi.org/10.1152/jappl.1999.86.5.1445
- [76]Passive force enhancement in striated muscleJournal of Applied Physiology 126:1782–1789https://doi.org/10.1152/japplphysiol.00676.2018
- [77]CNS learns stable, accurate, and efficient movements using a simple algorithmJournal of neuroscience 28:11–11https://doi.org/10.1523/JNEUROSCI.3099-08.2008
- [78]Signals & Systems
- [79]The spectral analysis of time series
- [80]Power amplification increases with contraction velocity during stretch-shortening cycles of skinned muscle fibersFrontiers in Physiology 12https://doi.org/10.3389/fphys.2021.644981
- [81]Molecular motors: Structural adaptations to cellular functionsNature 389:561–567https://doi.org/10.1038/39247
- [82]Titin-based contribution to shortening velocity of rabbit skeletal myofibrilsThe Journal of Physiology 540:177–188https://doi.org/10.1113/jphysiol.2001.013154
- [83]Mechanism of force enhancement during and after lengthening of active muscle: A temperature dependence studyJournal of muscle research and cell motility 33:313–325https://doi.org/10.1007/s10974-012-9307-8
- [84]Length dependence of changes in sarcoplasmic calcium concentration and myofibrillar calcium sensitivity in striated muscle fibresJournal of Muscle Research & Cell Motility 5:243–272https://doi.org/10.1007/BF00713107
- [85]The maximum shortening velocity of muscle should be scaled with activationJournal of Applied Physiology 86:1025–1031https://doi.org/10.1152/jappl.1999.86.3.1025
- [86]Single myosin molecule mechanics: Piconewton forces and nanometre stepsNature 368:113–119https://doi.org/10.1038/368113a0
- [87]Folding-unfolding transitions in single titin molecules characterized with laser tweezersScience 276:1112–1116https://doi.org/10.1126/science.276.5315.1112
- [88]Altered mechanical properties of titin immunoglobulin domain 27 in the presence of calciumEuropean Biophysics Journal 42:301–307https://doi.org/10.1007/s00249-012-0875-8
- [89]The role of sarcomere length non-uniformities in residual force enhancement of skeletal muscle myofibrilsRoyal Society Open Science 3https://doi.org/10.1098/rsos.150657
- [90]Biophysical Chemistry: Part I: The Conformation of Biological Macromolecules (Biophysical Chemistry)
- [91]The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodogramsIEEE Transactions on Audio and Electroacoustics 15:70–73https://doi.org/10.1109/TAU.1967.1161901
Article and author information
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:
- Reviewed Preprint version 3:
Copyright
© 2023, Millard et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
- views
- 232
- downloads
- 4
- citations
- 0
Views, downloads and citations are aggregated across all versions of this paper published by eLife.