The stiffness and damping of muscle are properties of fundamental importance for motor control, and the accurate simulation of muscle force. The CNS exploits the activation-dependent stiffness and damping (impedance) of muscle when learning new movements (Franklin et al., 2003), and when moving in unstable (Trumbower et al., 2009) or noisy environments (Selen et al., 2009). 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 (Burdet et al., 2001). With time and repetition, the force field becomes learned and co-contraction is reduced (Franklin et al., 2003).

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., 1994 were able to show that muscle behaves like a linear-time-invariant (LTI) system in response to small perturbations (see Appendix 9, Note 1): 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–90 Hz). When active muscle is stretched appreciably, titin can develop enormous forces (Herzog and Leonard, 2002; Leonard et al., 2010), which may prevent further lengthening and injury. The stiffness that best captures the response of muscle to the small perturbations of Kirsch et al., 1994 is far greater than the stiffness that best captures the response of muscle to large perturbations (Herzog and Leonard, 2002; Leonard et al., 2010). 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 (Kirsch et al., 1994) and large perturbations (Herzog and Leonard, 2002; Leonard et al., 2010) while also retaining the capability to reproduce the experiments of Hill, 1938 and Gordon et al., 1966. Unfortunately, this means that simulation studies that depend on an accurate replication of the perturbation response may reach conclusions well justified in simulation but not in reality. In this work, we focus on formulating a mechanistic muscle model (see Appendix 9, Note 2), 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 the famous force-velocity relationship of Hill, 1938, mechanistic Huxley (Huxley, 1957; Huxley and Simmons, 1971; Kosta et al., 2022) models in which individual elastic crossbridges are incorporated, and linearized muscle models (Hogan, 1985; Mussa-Ivaldi et al., 1985) which are accurate for small changes in muscle length. Kirsch et al., 1994 demonstrated that, for small perturbations, the force response of muscle is well represented by an activation-dependent spring and damper that are connected in parallel. Neither Hill nor Huxley models are likely to replicate the experiments of Kirsch et al., 1994 because a Hill muscle model (Zajac, 1989; Millard et al., 2013) does not contain any active spring elements; while a Huxley model lacks an active damping element. Although linearized muscle models can replicate the experiment of Kirsch et al., 1994, these models are only accurate for small changes in length and cannot replicate the nonlinear force-velocity relation of Hill, 1938, nor the nonlinear force-length relation of Gordon et al., 1966. 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., 1998 modeled the velocity dependence of muscle using a rheological element (see Appendix 9, Note 3) and an elastic rack rather than embedding the force-velocity relationship in equations directly, as is done in a typical Hill model (Zajac, 1989; Millard et al., 2013). This modification allows the model of Forcinito et al., 1998 to more faithfully replicate the force development of active muscle, as compared to a Hill-type model, during ramp length changes of ≈10% (see Appendix 9, Note 4) of the optimal CE length, and across velocities of 4–11% of the maximum contraction velocity (see Appendix 9, Note 5). Tamura and Saito, 2002 extended the work of Forcinito et al., 1998 by formulating a rheological muscle model with two Maxwell elements (a spring-damper in series) where one develops force quickly (high stiffness) and the other develops force slowly (low stiffness). By carefully selecting the dynamics that drive the two elements, the model of Tamura and Saito, 2002 replicated the force-length-velocity relations (Gordon et al., 1966; Hill, 1938) as well as qualitatively reproducing the force and stiffness profiles (Tamura et al., 2005) of force-enhancement and force-depression as measured by Sugi and Tsuchiya, 1988. Haeufle et al., 2014 made use of a serial-parallel network of spring-dampers to allow their model to reproduce the force-velocity relationship of Hill, 1938 mechanistically rather than embedding the experimental curve directly in their model. Günther et al., 2018 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 the experiment of Kirsch et al., 1994 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., 1994 specifically showed (see Figure 3 of Kirsch et al., 1994) that a serial connection of a spring-damper fails to reproduce the phase shift between force and length present in their experimental data.

Titin (Maruyama, 1976; Wang et al., 1979) has been more recently investigated to explain how lengthened muscle can develop active force when lengthened both within, and beyond, actin-myosin overlap (Leonard et al., 2010). 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 (Maruyama et al., 1985). 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 (Trombitás et al., 1998b) since the distal immunoglobulin (IgD) segments bind strongly to the thick filament (Houmeida et al., 1995). Titin has proven to be a complex filament, varying in composition and geometry between different muscle types (Tomalka et al., 2019; Boldt et al., 2020), widely between species (Lindstedt and Nishikawa, 2017), and can apply activation dependent forces to actin (Kellermayer and Granzier, 1996). 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 (Astier et al., 1998; Niederländer et al., 2004), between titin’s PEVK region and actin (Bianco et al., 2007; Nagy et al., 2004), between titin’s N2A region and actin (Dutta et al., 2018), and between the PEVK-IgD regions of titin and myosin (DuVall et al., 2017). 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 (Rode et al., 2009; Nishikawa et al., 2012; Schappacher-Tilp et al., 2015; Tahir et al., 2018; Heidlauf et al., 2016; Heidlauf et al., 2017), and more recently with myosin (DuVall et al., 2016).

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 ${\displaystyle 1.62{\mathcal{l}}_{\mathrm{o}}^{\mathrm{M}}}$ or less, the stiffness of the actin-myosin load path with a single attached crossbridge (0.22–1.15pN/nm) equals or exceeds the stiffness of 6 passive titin filaments (0.0348–0.173pN/nm), and our estimated stiffness of 6 active titin filaments (0.0696–0.346pN/nm, see Appendix 1 for further details). When fully activated, the stiffness of the actin-myosin load path (4.05–18.4pN/nm) far exceeds that of both the passive titin (0.0348–0.173pN/nm), and our estimated active titin (${\displaystyle 0.0696-0.346\text{pN/nm}}$) load paths. Since titin-focused models have not made any changes to the modeled myosin-actin interaction beyond a Hill (Zajac, 1989; Millard et al., 2013) or Huxley (Huxley, 1957; Huxley and Simmons, 1971) model, it is unlikely that these models would be able to replicate the experiments of Kirsch et al., 1994.

Although most motor control simulations (Kearney et al., 1997; Perreault et al., 2004; Schouten et al., 2008; Trumbower et al., 2009; Mitrovic et al., 2010) make use of the canonical linearized muscle model, phenomenological muscle models have also been used and modified to include stiffness. Sartori et al., 2015 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 (Kirsch et al., 1994), and there is no evidence for negative stiffness on the descending limb of the force-length curve (Herzog and Leonard, 2002). 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 (Heidlauf et al., 2017), 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., 2017 and De Groote et al., 2018 modeled short-range-stiffness using a stiff spring in parallel with the active force element of a Hill-type muscle model. While the approach developed in De Groote et al., 2017 and De Groote et al., 2018 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 (Kirsch et al., 1994, see Figure 9A); the damping of the canonical Hill-model has been left unchanged and likely differs substantially from biological muscle (Kirsch et al., 1994).

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 (Kirsch et al., 1994), while over longer time-scales the force-velocity (Hill, 1938) and force-length (Gordon et al., 1966) 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 the models of Rode et al., 2009 and (Schappacher-Tilp et al., 2015), features an activation-dependent (see Appendix 9, Note 6) 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 (van Soest et al., 2019), 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 bandwidth by simulating the stochastic perturbation experiments of Kirsch et al., 1994. The active-lengthening experiments of Herzog and Leonard, 2002 are used to evaluate the response of the model when it is actively lengthened within actin-myosin overlap. Next, we use the active-lengthening experiments of Leonard et al., 2010 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, 1938 and force-length experiments of Gordon et al., 1966. Since Hill-type models are so commonly used, we also replicate all of the simulated experiments using the Hill-type muscle model of Millard et al., 2013 to make the differences between these two types of models clear.

We begin by treating whole muscle as a scaled half-sarcomere that is pennated at an angle ${\displaystyle \alpha }$ with respect to a tendon (Figure 1A). The assumption that mechanical properties scale with size is commonly used when modeling muscle (Zajac, 1989) and makes it possible to model vastly different musculotendon units (MTUs) by simply changing the architectural and contraction properties: the maximum isometric force ${\displaystyle f_{\text{o}}^{\text{M}}}$, the optimal CE length ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ (at which the CE develops ${\displaystyle f_{\text{o}}^{\text{M}}}$), the pennation angle ${\displaystyle {\alpha }_{\text{o}}}$ of the CE (at a length of ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$) with respect to the tendon, the maximum shortening velocity ${\displaystyle v_{max}^{\text{M}}}$ of the CE, and the slack length of the tendon ${\displaystyle {\mathcal{l}}_{\text{s}}^{\text{T}}}$. Many properties of sarcomeres scale with ${\displaystyle f_{\text{o}}^{\text{M}}}$ and ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$: ${\displaystyle f_{\text{o}}^{\text{M}}}$ scales with physiological cross-sectional area (Maganaris et al., 2001), the force-length property scales with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ (Winters et al., 2011), the maximum normalized shortening velocity of different CE types scales with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ across animals great and small (Rome et al., 1990), and titin’s passive-force-length properties scale from single molecules to myofibrils (Herzog et al., 2012; Prado et al., 2005).

Figure 1

Download asset Open asset

The proposed model has several additional properties that we assume scale with ${\displaystyle f_{\text{o}}^{\text{M}}}$ and inversely with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$: the maximum active isometric stiffness ${\displaystyle k_{\text{o}}^{\text{X}}}$ and damping ${\displaystyle {\beta }_{\text{o}}^{\text{X}}}$, the passive forces due to the extracellular matrix (ECM), and passive forces due to titin. As crossbridge stiffness is well studied (Kaya and Higuchi, 2013), we assume that muscle stiffness due to crossbridges scales such that

where ${\displaystyle {\tilde{k}}_{\text{o}}^{\text{X}}}$ is the maximum normalized stiffness. This scaling is just what would be expected when many crossbridges (Kaya and Higuchi, 2013) act in parallel across the cross-sectional area of the muscle, and act in series along the length of the muscle. Although the intrinsic damping properties of crossbridges are not well studied, we assume that the linear increase in damping with activation observed by Kirsch et al., 1994 is due to the intrinsic damping properties of individual crossbridges which will also scale linearly with ${\displaystyle f_{\text{o}}^{\text{M}}}$ and inversely with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$

where ${\displaystyle {\tilde{\beta }}_{\text{o}}^{\text{X}}}$ is the maximum normalized damping. For the remainder of the paper, we refer to the proposed model as the VEXAT model due to the viscoelastic (VE) crossbridge (X) and active-titin (AT) elements of the model.

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 ${\displaystyle {\mathcal{l}}^{\text{S}}}$ (Figure 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 (Figure 1B). In addition, since the stiffness of actin (Higuchi et al., 1995) and myosin filaments (Tajima et al., 1994) greatly exceeds that of crossbridges (Veigel et al., 1998), 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 (Figure 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 (Figure 1C) and an additional passive load-path (Figure 1D). As is typical (Zajac, 1989; Millard et al., 2013), we assume that the passive elasticity of these structures scale linearly with ${\displaystyle f_{\text{o}}^{\text{M}}}$ and inversely with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. Since the VEXAT model has two passive load paths (Figure 1D), we further assume that the proportion of the passive force due to the extra-cellular-matrix (ECM) and titin does not follow a scale-dependent pattern, but varies from muscle-to-muscle as observed by Prado et al., 2005.

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 (Dutta et al., 2018), there is also support for a titin-actin interaction occurring near titin’s PEVK region (Bianco et al., 2007; Nagy et al., 2004), and for a titin-myosin interaction near the PEVK-IgD region (DuVall et al., 2017). 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 (Bianco et al., 2007; Nagy et al., 2004) in the PEVK segment; and second, there is evidence supporting an interaction at the proximal end of the PEVK segment (Dutta et al., 2018, N2A-actin interaction). 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, 1996 and Dutta et al., 2018: 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 ${\displaystyle f_{\text{o}}^{\text{M}}}$ and inversely with ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$.

After lumping all of the crossbridges and titin filaments together, we are left with a rigid-tendon MTU model that has two generalized positions

and an elastic-tendon MTU model that has three generalized positions

Given these generalized positions, the path length ${\displaystyle {\mathcal{l}}^{\text{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 (Figure 1A) to prevent the model from reaching the numerical singularity that exists as ${\displaystyle {\tilde{\mathcal{l}}}^{\text{M}}}$ approaches ${\displaystyle {\tilde{\mathcal{l}}}_{min}^{\text{M}}}$, the length at which ${\displaystyle \alpha \to {90}^{\circ }}$ in Equation 6 and Equation 7. The tendon length

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 ${\displaystyle {\mathcal{l}}^{\text{S}}}$ and the length of the myosin filament ${\displaystyle {\text{L}}^{\text{M}}}$. The length of ${\displaystyle {\mathcal{l}}^{\text{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 (${\displaystyle {\mathcal{l}}^{\text{1}}+{\text{L}}^{\text{T12}}}$) and the length of the IgD segment that is bound to myosin (${\displaystyle {\text{L}}^{\text{IgD}}}$). Finally, the length of the extra-cellular-matrix ${\displaystyle {\mathcal{l}}^{\text{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: the stochastic perturbations of Kirsch et al., 1994 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 ${\displaystyle {\mathcal{l}}^{\text{S}}}$, and velocity ${\displaystyle v^{\text{S}}}$. When the activation state ${\displaystyle a}$ and the titin-actin interaction model are included, the resulting rigid-tendon model 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 (${\displaystyle a}$), while an elastic-tendon model has two states (${\displaystyle a}$ and ${\displaystyle {\mathcal{l}}^{\text{M}}}$) (Millard et al., 2013).

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 ${\displaystyle C^2}$ continuous (see Appendix 9, Note 7) 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 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$, forces by the maximum active isometric tension ${\displaystyle f_{\text{o}}^{\text{M}}}$, stiffness and damping by ${\displaystyle f_{\text{o}}^{\text{M}}/{\mathcal{l}}_{\text{o}}^{\text{M}}}$. Velocities used as input to the force-velocity relation ${\displaystyle {\mathbf{f}}^{\text{V}}}$ are further normalized by ${\displaystyle v_{\text{max}}^{\text{M}}}$ and annotated using a hat: ${\displaystyle {\hat{v}}^{\text{M}}=v^{\text{M}}/v_{\text{max}}^{\text{M}}}$. Tendon lengths and velocities are normalized by ${\displaystyle {\mathcal{l}}_{\text{s}}^{\text{T}}}$ tendon slack length, while forces are normalized by ${\displaystyle f_{\mathrm{o}}^{\mathrm{M}}}$.

To evaluate the state derivative of the model, we require equations for ${\displaystyle \dot{a}}$, ${\displaystyle {\dot{v}}^{\text{S}}}$, ${\displaystyle v^1}$, and ${\displaystyle v^{\text{M}}}$ if the tendon is elastic. For ${\displaystyle \dot{a}}$ we use of the first order activation dynamics model described in Millard et al., 2013 (see Appendix 9, Note 8) which uses a lumped first-order ordinary-differential-equation (ODE) to describe how a fused tetanus electrical excitation leads to force development in an isometric muscle. We formulated the equation for ${\displaystyle {\dot{v}}^{\text{S}}}$ with the intention of having the model behave like a spring-damper over small time-scales, but to converge to the tension developed by a Hill-type model

over longer time-scales, where ${\displaystyle {\mathbf{f}}^{\text{L}}(\cdot )}$ is the active-force-length curve (Figure 2A), ${\displaystyle {\mathbf{f}}^{\text{PE}}(\cdot )}$ is the passive-force-length curve (Figure 2A), and ${\displaystyle {\mathbf{f}}^{\text{V}}(\cdot )}$ is the force-velocity curve (Figure 2B).

Figure 2

Download asset Open asset

The normalized tension developed by the VEXAT model

differs from that of a Hill model, Equation 14, because it has no explicit dependency on ${\displaystyle {\tilde{v}}^{\text{M}}}$, includes four passive terms, and a lumped viscoelastic crossbridge element. The four passive terms come from the ECM element ${\displaystyle {\mathbf{f}}^{\text{ECM}}({\tilde{\mathcal{l}}}^{\text{ECM}})}$ (Figure 3A), the PEVK-IgD element ${\displaystyle {\mathbf{f}}^2({\tilde{\mathcal{l}}}^2)}$ (Figure 3A and B), the compressive term ${\displaystyle {\mathbf{f}}^{\text{KE}}({\tilde{\mathcal{l}}}^{\text{M}})}$ (prevents ${\displaystyle {\tilde{\mathcal{l}}}^{\text{M}}\cos \alpha }$ from reaching a length of 0), and a numerical damping term ${\displaystyle {\tilde{\beta }}^{ϵ}{\tilde{v}}^{\text{M}}}$ (where ${\displaystyle {\tilde{\beta }}^{ϵ}}$ is small). The active force developed by the XE’s lumped crossbridge ${\displaystyle {\tilde{k}}_{\text{o}}^{\text{X}}{\tilde{\mathcal{l}}}^{\text{X}}+{\tilde{\beta }}_{\text{o}}^{\text{X}}{\tilde{v}}^{\text{X}}}$ is scaled by the fraction of the XE that is active and attached, ${\displaystyle a{\mathbf{f}}^{\text{L}}({\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}})}$, where ${\displaystyle {\mathbf{f}}^{\text{L}}(\cdot )}$ is the active-force-length relation (Figure 2A). We evaluate ${\displaystyle {\mathbf{f}}^{\text{L}}}$ using ${\displaystyle {\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}}}$ in Equation 15, rather than ${\displaystyle {\tilde{\mathcal{l}}}^{\text{M}}}$ as in Equation 14, since actin-myosin overlap is independent of crossbridge strain. With ${\displaystyle {\tilde{f}}^{\text{M}}}$ derived, we can proceed to model the acceleration of the CE, ${\displaystyle {\dot{v}}^{\text{S}}}$, so that it is driven over time by the force imbalance between the XE’s active tension and that of a Hill model.

Figure 3

Download asset Open asset

We set the first term of ${\displaystyle {\dot{\tilde{v}}}^{\text{S}}}$ so that, over time, the CE is driven to develop the same active tension as the Hill-type model of Millard et al., 2013 (terms highlighted in blue)

where ${\displaystyle {\tau }^{\text{S}}}$ is a time constant and ${\displaystyle {\mathbf{f}}^{\text{V}}({\hat{v}}^{\text{S}})}$ is the force-velocity curve (Figure 2B). The rate of adaptation of the model’s tension, to the embedded Hill model, is set by the time constant ${\displaystyle {\tau }^{\text{S}}}$: as ${\displaystyle {\tau }^{\text{S}}}$ is decreased the VEXAT model converges more rapidly to a Hill-type model; as ${\displaystyle {\tau }^{\text{S}}}$ is increased the active force produced by the model will look more like a spring-damper. Our preliminary simulations indicate that there is a trade-off to choosing ${\displaystyle {\tau }^{\text{S}}}$: when ${\displaystyle {\tau }^{\text{S}}}$ is large the model will not shorten rapidly enough to replicate Hill’s experiments, while if ${\displaystyle {\tau }^{\text{S}}}$ is small the low-frequency response of the model is compromised when the experiments of Kirsch et al., 1994 are simulated.

The remaining two terms, ${\displaystyle \text{D}{\tilde{v}}^{\text{S}}}$ and ${\displaystyle {\text{e}}^{-(a/a_{\text{L}}{)}^2}({\text{G}}_{\text{L}}{\tilde{\mathcal{l}}}^{\text{X}}+{\text{G}}_{\text{V}}{\tilde{v}}^{\text{X}})}$, have been included for numerical reasons specific to this model formulation rather than muscle physiology. We include a term that damps the rate of actin-myosin translation, ${\displaystyle \text{D}{\tilde{v}}^{\text{S}}}$, to prevent this second-order system from unrealistically oscillating (see Appendix 9, Note 9). The final term ${\displaystyle {\text{e}}^{-(a/a_{\text{L}}{)}^2}({\text{G}}_{\text{L}}{\tilde{v}}^{\text{X}}+{\text{G}}_{\text{V}}{\tilde{v}}^{\text{X}})}$, where ${\displaystyle {\text{G}}_{\text{L}}}$ and ${\displaystyle {\text{G}}_{\text{V}}}$ are scalar gains, and ${\displaystyle a_{\text{L}}}$ is a low-activation threshold (${\displaystyle a_{\text{L}}}$ is ${\displaystyle 0.05}$ in this work). This final term has been included as a consequence of the generalized positions we have chosen. When the CE is nearly deactivated (as ${\displaystyle a}$ approaches ${\displaystyle a_{\text{L}}}$), this term forces ${\displaystyle {\tilde{\mathcal{l}}}^{\text{S}}}$ and ${\displaystyle {\tilde{v}}^{\text{S}}}$ to shadow the location and velocity of the XE attachment point. This ensures that if the XE is suddenly activated, that it attaches with little strain. We had to include this term because we made ${\displaystyle {\mathcal{l}}^{\text{S}}}$ a state variable, rather than ${\displaystyle {\mathcal{l}}^{\text{X}}}$. We chose ${\displaystyle {\mathcal{l}}^{\text{S}}}$ as a state variable, rather than ${\displaystyle {\mathcal{l}}^{\text{X}}}$, so that the states are more equally scaled for numerical integration.

The passive force developed by the CE in Equation 15 is the sum of the elastic forces (Figure 3A) developed by the force-length curves of titin (${\displaystyle {\mathbf{f}}^{\text{1}}({\tilde{\mathcal{l}}}^{\text{1}})}$ and ${\displaystyle {\mathbf{f}}^{\text{2}}({\tilde{\mathcal{l}}}^{\text{2}})}$) and the ECM (${\displaystyle {\mathbf{f}}^{\text{ECM}}({\tilde{\mathcal{l}}}^{\text{ECM}})}$). We model titin’s elasticity as being due to two serially connected elastic segments: the first elastic segment ${\displaystyle {\mathbf{f}}^1({\tilde{\mathcal{l}}}^1)}$ is formed by lumping together the IgP segment and a fraction ${\displaystyle \text{Q}}$ of the PEVK segment, while the second elastic segment ${\displaystyle {\mathbf{f}}^2({\tilde{\mathcal{l}}}^2)}$ is formed by lumping together the remaining ${\displaystyle (1-\text{Q})}$ of the PEVK segment with the free IgD section. Our preliminary simulations of the active lengthening experiment of Herzog and Leonard, 2002 indicate that a ${\displaystyle \text{Q}}$ value of 0.5, positioning the PEVK-actin attachment point that is near the middle of the PEVK segment, allows the model to develop sufficient tension when actively lengthened. The large section of the IgD segment that is bound to myosin is treated as rigid.

The curves that form ${\displaystyle {\mathbf{f}}^{\text{ECM}}({\tilde{\mathcal{l}}}^{\text{ECM}})}$, ${\displaystyle {\mathbf{f}}^1({\tilde{\mathcal{l}}}^1)}$, and ${\displaystyle {\mathbf{f}}^2({\tilde{\mathcal{l}}}^2)}$ have been carefully constructed to satisfy three experimental observations: that the total passive force-length curve of titin and the ECM match the observed passive force-length curve of the muscle (Figure 2A and Figure 3A) as in the experiments of Prado et al., 2005; that the proportion of the passive force developed by titin and the ECM (Figure 3A) is within experimental observations of Prado et al., 2005; and that the Ig domains and PEVK residues show the same relative elongation (Figure 3C) as observed by Trombitás et al., 1998a. Even though the measurements of Trombitás et al., 1998b come from human soleus titin, we can construct the force-length curves of other titin isoforms if the number of proximal Ig domains, PEVK residues, and distal Ig domains are known (see Appendix 2.4). Since the passive–force-length relation and the results of Trombitás et al., 1998b are at modest lengths, we consider two different extensions to the force-length relation to simulate extreme lengths: first, a simple linear extrapolation; second, we extend the force-length relation of each of titin’s segments to follow a worm-like-chain (WLC) model similar to Trombitás et al., 1998b (see Appendix 2.4 for details on the WLC model). With titin’s passive force-length relations defined, we can next consider what happens to titin in active muscle.

When active muscle is lengthened, it produces an enhanced force that persists long after the lengthening has ceased called residual force enhancement (RFE) (Herzog and Leonard, 2002). For the purposes of the VEXAT model, we assume that RFE is produced by titin. Experiments have demonstrated RFE on both the ascending limb (Peterson et al., 2004) and descending limb of the force-length (Herzog and Leonard, 2002) relation. The amount of RFE depends both on the final length of the stretch (Hisey et al., 2009) 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, ${\displaystyle {\mathcal{l}}^{\text{2}}}$, can contribute to RFE when the titin-actin bond is formed and CE is lengthened beyond ${\displaystyle {\tilde{\mathcal{l}}}_s^{\text{M}}}$, the shortest CE length at which the PEVK-actin bond can form. In this work, we set ${\displaystyle {\tilde{\mathcal{l}}}_s^{\mathrm{M}}}$ to be equal to the slack length of the CE’s force-length relation ${\displaystyle {\tilde{\mathcal{l}}}_{\text{s}}^{\mathrm{P}\mathrm{E}}}$ (see Table 1E and H). To incorporate the asymmetric length dependence of RFE observed by Hisey et al., 2009, we introduce a smooth step-up function

that transitions from zero to one as ${\displaystyle {\tilde{\mathcal{l}}}^{\text{M}}}$ extends beyond ${\displaystyle {\tilde{\mathcal{l}}}_s^{\text{M}}}$, where the sharpness of the transition is controlled by ${\displaystyle \text{R}}$. Similar to the experimental work of Hisey et al., 2009, active lengthening on the ascending limb will produce similar amounts of RFE since ${\displaystyle {\tilde{\mathcal{l}}}_s^{\text{M}}<{\mathcal{l}}_{\text{o}}^{\text{M}}}$ and the titin-actin bond is prevented from forming below ${\displaystyle {\tilde{\mathcal{l}}}_s^{\text{M}}}$. In contrast, the amount of RFE on the descending limb increases with increasing stretch magnitudes since the titin-actin bond is able to form across the entire descending limb.

At very long CE lengths, the modeled titin-actin bond can literally slip off of the end of the actin filament (Figure 1A) when the distance between the Z-line and the bond, ${\displaystyle {\tilde{\mathcal{l}}}^{\text{1}}+{\tilde{\text{L}}}^{\text{T12}}}$, exceeds the length of the actin filament, ${\displaystyle {\tilde{\text{L}}}^{\text{A}}}$. To break the titin-actin bond at long CE lengths, we introduce a smooth step-down function

The step-down function ${\displaystyle u^{\text{L}}}$ transitions from one to zero when the titin-actin bond (${\displaystyle {\tilde{\mathcal{l}}}^{\text{1}}+{\tilde{\text{L}}}^{\text{T12}}}$) reaches ${\displaystyle {\tilde{\text{L}}}^{\text{A}}}$, the end of the actin filament.

The strength of the titin-actin bond also appears to vary nonlinearly with activation. Fukutani and Herzog, 2018 observed that the absolute RFE magnitude produced by actively lengthened fibers is similar between normal and reduced contractile force states. Since the experiments of Fukutani and Herzog, 2018 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, 2018 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 ${\displaystyle {\text{A}}_{\circ }}$ 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 Ca²⁺ concentration (Millar and Homsher, 1990), this is a mathematical convenience. The work of Leonard et al., 2010 makes it clear that both Ca²⁺ and crossbridge cycling are needed to allow titin to develop enhanced forces during active lengthening: no enhanced forces are observed in the presence of Ca²⁺ 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 ${\displaystyle u^{\text{A}}{u}^{\text{S}}{u}^{\text{L}}{\beta }_{\text{A}}^{\text{PEVK}}}$, where ${\displaystyle {\beta }_{\text{A}}^{\text{PEVK}}}$ is the maximum damping coefficient.

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 (${\displaystyle u^{\text{A}}{u}^{\text{S}}{u}^{\text{L}}{\beta }_{\text{A}}^{\text{PEVK}}}$) and a small amount of passive damping (${\displaystyle {\beta }_{\text{P}}^{\text{PEVK}}}$). Since Equation 20 is linear in ${\displaystyle {\tilde{v}}^{\text{1}}}$, we can solve directly for it

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 Equations 6–11 and the generalized positions, evaluating the generalized velocities of a rigid-tendon and elastic-tendon model differ substantially. The CE velocity ${\displaystyle v^{\text{M}}}$ and pennation angular velocity ${\displaystyle \dot{\alpha }}$ of a rigid-tendon model can be evaluated directly given the path length, velocity, and the time derivatives of Equation 6 and Equation 8. After ${\displaystyle v^{\text{1}}}$ is evaluated using Equation 21, the velocities of the remaining segments can be evaluated using the time derivatives of Equations 9–11.

Evaluating the CE rate of lengthening, ${\displaystyle v^{\text{M}}}$, for an elastic-tendon muscle model is more involved. As is typical of lumped parameter muscle models (Zajac, 1989; Thelen, 2003; Millard et al., 2013), here we assume that difference in tension, ${\displaystyle {\tilde{f}}^{ϵ}}$, between the CE and the tendon

is negligible (see Appendix 9, Note 10). 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., 1994 particularly at high frequencies. Kirsch et al., 1994 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 (see Appendix 9, Note 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 ${\displaystyle U{\hat{\mathbf{k}}}^{\text{T}}({\tilde{\mathcal{l}}}^{\text{T}})}$, is a linear scaling of the normalized tendon stiffness ${\displaystyle {\hat{\mathit{k}}}^{\text{T}}}$ by ${\displaystyle U}$, a constant scaling coefficient. We have chosen this specific damping model because it fits the data of Netti et al., 1996 and captures the structural coupling between tendon stiffness and damping (see Appendix 2.2 and Figure 1 for further details).

Now that all the terms in Equation 22 have been explicitly defined, we can use Equation 22 to solve for ${\displaystyle v^{\text{M}}}$. Equation 22 becomes linear in ${\displaystyle v^{\text{M}}}$ after substituting the force models described in Equation 23 and Equation 15, and the kinematic model described in Equation 8, Equation 9 and Equation 11 (along with the time derivatives of Equations 8–11). After some simplification we arrive at

allowing us to evaluate the final state derivative in ${\displaystyle \underset{\_}{\dot{x}}}$. During simulation the denominator of ${\displaystyle {\tilde{v}}^{\text{M}}}$ will always be finite since ${\displaystyle {\tilde{\beta }}^{ϵ}>0}$, and ${\displaystyle \alpha <{90}^{\circ }}$ due to the compressive element. The evaluation of ${\displaystyle \underset{\_}{\dot{x}}}$ in the VEXAT model is free of numerical singularities, giving it an advantage over a conventional Hill-type muscle models (Millard et al., 2013). In addition, the VEXAT’s ${\displaystyle \underset{\_}{\dot{x}}}$ does not require iteration to numerically solve a root, giving it an advantage over a singularity-free formulation of the Hill model (Millard et al., 2013). As with previous models, initializing the model’s state is not trivial and required the derivation of a model-specific method (see Appendix 3 for details).

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 (1–3.8% ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$) stochastic perturbation experiment of Kirsch et al., 1994 demonstrates that the impedance of muscle is well described by a stiff spring in parallel with a damper, and that the spring-damper coefficients vary linearly with active force. The active impedance of muscle is such a fundamental part of motor learning that the amount of impedance, as indicated by co-contraction, is used to define how much learning has actually taken place (Franklin et al., 2003; Franklin et al., 2008): co-contraction is high during initial learning, but decreases over time as a task becomes familiar. The active lengthening experiment of Herzog and Leonard, 2002 shows that modestly stretched (7–21% ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$) biological muscle has positive stiffness even on the descending limb of the active force-length curve (${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}>1}$). In contrast, a conventional Hill model (Zajac, 1989; Millard et al., 2013) can have negative stiffness on the descending limb of the active-force-length curve, a property that is both mechanically unstable and unrealistic. The final active lengthening experiment of Leonard et al., 2010 unequivocally demonstrates that the CE continues to develop active forces during extreme lengthening (329% ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$) which exceeds actin-myosin overlap. Active force development beyond actin-myosin overlap is made possible by titin, and its activation-dependent interaction with actin (Leonard et al., 2010). The biological benchmark simulations conclude with a replication of the force-velocity experiments of Hill, 1938 and the force-length experiments of Gordon et al., 1966.

The VEXAT model requires the architectural muscle parameters (${\displaystyle f_{\text{o}}^{\text{M}}}$, ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$, ${\displaystyle {\alpha }_{\text{o}}}$, ${\displaystyle v_{max}^{\text{M}}}$, and ${\displaystyle {\mathcal{l}}_{\text{s}}^{\text{T}}}$) needed by a conventional Hill-type muscle model as well as additional parameters. The additional parameters are needed for these component models: the compressive element (Equation 15 and Equation 24), the lumped viscoelastic XE (Equation 1 and Equation 2), XE-actin dynamics (Equation 16), the two-segment active titin model (Figure 3), titin-actin dynamics (Equation 21), and the tendon damping model (Equation 23). Fortunately, there is enough experimental data in the literature that values can be found, fitted, or estimated directly for our simulations of experiments on cat soleus (Table 1F.-I.), and rabbit psoas fibrils (see Appendix 2 and Appendix 8 for rabbit psoas fibril model parameters). The parameter values we have established for the cat soleus (Table 1F.-I.) can serve as initial values when modeling other mammalian MTU’s because these parameters have been normalized (by ${\displaystyle f_{\text{o}}^{\text{M}}}$, ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$, and ${\displaystyle {\mathcal{l}}_{\text{s}}^{\text{T}}}$ where appropriate) and will scale appropriately given the architectural properties of a different MTU. By making use of these default values, the VEXAT model can be made to represent another MTU using exactly the same number of parameters as a Hill-type muscle model (Table 1A.-E.).

Stochastic length perturbation experiments

In the in-situ experiment of Kirsch et al., 1994, 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., 1994 applied stochastic length perturbations (Figure 4A, blue line) to elicit force responses from the muscle (in this case a spring-damper Figure 4A, black line) across a broad range of frequencies (4–90 Hz) and across a range of small length perturbations (1–3.8% ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$). The resulting time-domain signals can be quite complicated (Figure 4A) but contain rich measurements of how muscle transforms changes in length into changes in force.

Figure 4

Download asset Open asset

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 (Oppenheim et al., 1996; Koopmans, 1995) can be applied to extract a relationship between length ${\displaystyle x(t)}$ and force ${\displaystyle y(t)}$. We will give a brief overview of system identification methods here to make methods and results clearer. First, the time-domain signals (${\displaystyle x(t)}$ and ${\displaystyle y(t)}$) are transformed into an equivalent representation in the frequency-domain (${\displaystyle X(s)}$ and ${\displaystyle Y(s)}$) as a sum of scaled and shifted sine curves (Figure 4B and C) using a Fourier transform (Oppenheim et al., 1996). In the frequency domain, we identify an LTI system of best fit ${\displaystyle H(s)}$ that describes how muscle transforms changes in length into changes in force such that ${\displaystyle Y(s)=H(s)X(s)}$. Next, we evaluate how ${\displaystyle H(s)}$ scales the magnitude (gain) and shifts the timing (phase) of a sinusoid in ${\displaystyle X(s)}$ into a sinusoid of the same frequency in ${\displaystyle Y(s)}$ (Figure 4D). This process is repeated across all frequency-matched pairs of input and output sinusoids to build a function of how muscle scales (Figure 4E) and shifts (Figure 4F) input length sinusoids into output force sinusoids. The resulting transformation turns two complicated time-domain signals (Figure 4A) into a clear relationship in the frequency-domain that describes how muscle transforms length changes into force changes: a very slow (≈0 Hz) length change will result in an output force that is scaled by 4.5 and is in phase (Figure 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° (Figure 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 Figure 4A, 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 ${\displaystyle H(s)}$ and the corresponding gain and phase profiles is more elaborate than we have described (see Appendix 4 for details).

Kirsch et al., 1994 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., 1994 used the coherence squared, ${\displaystyle (C_{xy}{)}^2}$, between the input and output to evaluate the degree of linearity: ${\displaystyle Y(s)}$ is a linear transformation of ${\displaystyle X(s)}$ at frequencies in which ${\displaystyle (C_{xy}{)}^2}$ is near one, while ${\displaystyle Y(s)}$ cannot be described as a linear function of ${\displaystyle X(s)}$ at frequencies in which ${\displaystyle (C_{xy}{)}^2}$ approaches zero. By calculating ${\displaystyle (C_{xy}{)}^2}$ between the length perturbation and force waveforms, Kirsch et al., 1994 identified the bandwidth in which the muscle’s response is approximately linear. Kirsch et al., 1994 set the lower frequency of this band to 4 Hz, and Figure 3 of Kirsch et al., 1994 suggests that this corresponds to ${\displaystyle (C_{xy}{)}^2\ge 0.67}$ though the threshold for ${\displaystyle (C_{xy}{)}^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., 1994 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, Kirsch et al., 1994 evaluated the stiffness and damping coefficients that best fit the muscle’s frequency response. Finally, Kirsch et al., 1994 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

(25) ${\displaystyle VAF(f^{\text{KD}},f^{\text{EXP}})=\frac{{\sigma }^2(f^{\text{EXP}})-{\sigma }^2(f^{\text{KD}}-f^{\text{EXP}})}{{\sigma }^2(f^{\text{EXP}})}.}$

Astonishingly, Kirsch et al., 1994 found that a spring-damper of best fit has a VAF of between 78% and 99% (see Appendix 9, Note 12) when compared to the experimentally measured forces ${\displaystyle f^{\text{EXP}}}$. By repeating this experiment over a variety of stimulation levels (using both electrical stimulation and the crossed-extension reflex) Kirsch et al., 1994 showed that these stiffness and damping coefficients vary linearly with the active force developed by the muscle. Furthermore, Kirsch et al., 1994 repeated the experiment using perturbations that had a variety of lengths (0.4 mm, 0.8 mm, and 1.6 mm) and bandwidths (15 Hz, 35 Hz, and 90 Hz) 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., 1994 for details). Here, we simulate the experiment of Kirsch et al., 1994 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., 1994 studied.

To simulate the experiments of Kirsch et al., 1994 we begin by creating the 9 stochastic perturbation waveforms used in the experiment that vary in perturbation amplitude (0.4 mm, 0.8mm, and 1.6 mm) and bandwidth (0–15 Hz, 0–35 Hz, and 0–90 Hz) (see Appendix 9, Note 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 second-order Butterworth filter is applied to the waveform and finally the signal is scaled to the appropriate amplitude (Figure 5). The muscle model is then activated until it develops a constant tension at a length of ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. The musculotendon unit is then simulated as the length is varied using the previously constructed waveforms while activation is held constant. To see how impedance varies with active force, we repeated these simulations at ten evenly spaced tensions from 2.5N to 11.5N. Ninety simulations are required to evaluate the nine different perturbation waveforms at each of the ten tension levels. The time-domain length perturbations and force responses of the modeled muscles are used to evaluate the coherence squared of the signal, gain response, and phase responses of the models in the frequency-domain. Since the response of the models might be more nonlinear than biological muscle, we select a bandwidth that meets ${\displaystyle (C_{xy}{)}^2>0.67}$ but otherwise follows the bandwidths analyzed by Kirsch et al., 1994 (see Appendix 4 for details).

Figure 5

Download asset Open asset

When coupled with an elastic-tendon, the 15 Hz perturbations show that neither model can match the VAF of the analysis of Kirsch et al., 1994 (compare Figure 6A to G), while at 90 Hz the VEXAT model reaches a VAF of 89% (Figure 6D) which is within the range of 78–99% reported by Kirsch et al., 1994. In contrast, the Hill model’s VAF at 90 Hz remains low at 58% (Figure 6J). While the VEXAT model has a gain profile in the frequency-domain that closer to the data of Kirsch et al., 1994 than the Hill model (compare Figure 6B to H and E to K), both models have a greater phase shift than the data of Kirsch et al., 1994 at low frequencies (compare Figure 6C to I and F to L). The phase response of the VEXAT model to the 90 Hz perturbation (Figure 6F) shows the consequences of Equation 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 Equation 16: the value of ${\displaystyle {\tau }^{\text{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 15 Hz and 90 Hz profiles for the VEXAT model (listed on Figure 6A. and D.) and the Hill model (listed on Figure 6G. and J., respectively).

Figure 6

Download asset Open asset

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 Figure 6A–F to Figure 1A–F in Appendix 6), the response of the Hill model to the 15 Hz perturbation changes dramatically with the time-domain VAF rising from 55% to 85% (compare Figure 6G–L to Figure 1G–L in Appendix 6). 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 (Figure 1H in Appendix 6), 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 (Kirsch et al., 1994).

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 Figure 6A–F to Figure 1A–F in Appendix 6), 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 (Figure 1G–L in Appendix 6). The stiffness of a Hill model’s CE

(26) ${\displaystyle k^{\text{M}}=f_{\text{o}}^{\text{M}}\left(a\frac{\text{d}{\mathbf{f}}^{\text{L}}}{\text{d}{\mathcal{l}}^{\text{M}}}{\mathbf{f}}^{\text{V}}+\frac{\text{d}{\mathbf{f}}^{\text{PE}}}{\text{d}{\mathcal{l}}^{\text{M}}}\right)}$

is heavily influenced by the partial derivative of ${\displaystyle \frac{\text{d}{\mathbf{f}}^{\text{L}}}{\text{d}{\mathcal{l}}^{\text{M}}}}$ which has a region of negative stiffness. Although ${\displaystyle \frac{\text{d}{\mathbf{f}}^{\text{PE}}}{\text{d}{\mathcal{l}}^{\text{M}}}}$ is well approximated as being linear for small length changes, ${\displaystyle \frac{\text{d}{\mathbf{f}}^{\text{L}}}{\text{d}{\mathcal{l}}^{\text{M}}}}$ changes sign across ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. The damping of a Hill model’s CE

(27) ${\displaystyle {\beta }^{\text{M}}=f_{\text{o}}^{\text{M}}\left(a{\mathbf{f}}^{\text{L}}\frac{\text{d}{\mathbf{f}}^{\text{V}}}{\text{d}{v}^{\text{M}}}\right)}$

also suffers from high degrees of nonlinearity for small perturbations about ${\displaystyle v^{\text{M}}=0}$ since the slope of ${\displaystyle \frac{\text{d}{\mathbf{f}}^{\text{V}}}{\text{d}{v}^{\text{M}}}}$ is positive and large when shortening, and positive and small when lengthening (Figure 2B). While Equation 26 and Equation 27 are mathematically correct, the negative stiffness and wide ranging damping values predicted by these equations do not match experimental data (Kirsch et al., 1994). In contrast, the stiffness

(28) ${\displaystyle {\tilde{k}}^{\text{M}}=\frac{a{\tilde{k}}_{\text{o}}^{\text{X}}{\mathbf{f}}^{\text{L}}({\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}})}{2}+\frac{\text{d}{\mathbf{f}}^{\text{2}}({\tilde{\mathcal{l}}}^{\text{2}})}{\text{d}{\mathcal{l}}^{\text{2}}}\frac{1}{2}+\frac{\text{d}{\mathbf{f}}^{\text{ECM}}({\tilde{\mathcal{l}}}^{\text{ECM}})}{\text{d}{\mathcal{l}}^{\text{ECM}}}\frac{1}{2}}$

and damping

(29) ${\displaystyle {\tilde{\beta }}^{\text{M}}=a{\mathbf{f}}^{\text{L}}({\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}})\left({\tilde{\beta }}_{\text{o}}^{\text{X}}\frac{\text{d}{\tilde{v}}^{\text{X}}}{\text{d}{\tilde{v}}^{\text{M}}}\right)+{\tilde{\beta }}^{ϵ}}$

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 2.5 for derivation).

By repeating the stochastic perturbation experiments across a range of isometric forces, Kirsch et al., 1994 were able to show that the stiffness and damping of a muscle varies linearly with the active tension it develops (see Figure 12 of Kirsch et al., 1994). We have repeated our simulations of the experiments of Kirsch et al., 1994 at ten nominal forces (spaced evenly between 2.5 N and 11.5 N) and compared how the VEXAT model and the Hill model’s stiffness and damping coefficients (Figure 7) compare to Figure 12 of Kirsch et al., 1994. The stiffness and damping profile of the VEXAT model deviates a little from the data of Kirsch et al., 1994 because XE’s dynamics at 35 Hz are still influenced by the Hill model embedded in Equation 16 (see Appendix 2.5). Despite this, the VEXAT model develops similar stiffness and damping profile with either a viscoelastic-tendon (Figure 7A and B) or a rigid-tendon (Figure 7C and D). In contrast, when the Hill model is coupled with an elastic-tendon both its stiffness and damping are larger than what is reported by Kirsch et al., 1994 at the higher tensions (Figure 7A and B). This pattern changes when simulating a Hill model with a rigid-tendon: the model’s stiffness is slightly negative (Figure 7C), while the model’s final damping coefficient is nearly three times the value measured by Kirsch et al., 1994 (Figure 7D). Though a negative stiffness may seem surprising, Equation 26 shows a negative stiffness is possible at the nominal CE length of these simulations: just past ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ the slope of the active force-length curve is negative and the slope of the passive force-length curve is negligible. The tendon model also affects the VAF of the Hill model to a large degree: the elastic-tendon Hill model has a low VAF 30–51% (Figure 7A & B) while the rigid-tendon Hill model has a much higher VAF of 86%. Although the VAF of the rigid-tendon Hill model is acceptable, these forces are being generated in a completely different manner than those obtained from biological muscle, as the data of Kirsch et al., 1994 indicate (Figure 7C and D).

Figure 7

Download asset Open asset

When the VAF of the VEXAT and Hill model is evaluated across a range of nominal tensions (ten values from 2.5 to 11.5 N), frequencies (15 Hz, 35 Hz, and 90 Hz), amplitudes (0.4 mm, 0.8 mm, and 1.6 mm), 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., 1994 while the Hill model’s 28–95% VAF differs (Figure 8); and second, that there are systematic variations in VAF, stiffness, and damping across the different perturbation magnitudes and frequencies (see Appendix 5—table 1 and Appendix 5—table 2). Both models produce worse VAF values when coupled with an elastic-tendon (Figure 8A, B and C), although 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 (Figure 8A). While the VEXAT model’s lowest VAF occurs in response to the low frequency perturbations (Figure 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 (Figure 8C) while the elastic-tendon Hill model has its lowest VAF in response to the 0.4 mm 15 Hz perturbations (Figure 8A). It is unclear if biological muscle displays systematic shifts in VAF since Kirsch et al., 1994 did not report the VAF of each trial.

Figure 8

Download asset Open asset

Active lengthening on the descending limb

We now turn our attention to the active lengthening in-situ experiments of Herzog and Leonard, 2002. During these experiments, cat soleus muscles were actively lengthened by modest amounts (7–21% ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$) starting on the descending limb of the active-force-length curve (${\displaystyle {\mathcal{l}}^{\text{M}}/{\mathcal{l}}_{\text{o}}^{\text{M}}>1}$ in Figure 2A). This starting point was chosen specifically because the stiffness of a Hill model may actually change sign and become negative because of the influence of the active-force-length curve on ${\displaystyle k^{\text{M}}}$ as shown in Equation 26 as ${\displaystyle {\mathcal{l}}^{\text{M}}}$ extends beyond ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. The experiment of Herzog and Leonard, 2002 is important for showing that biological muscle does not exhibit negative stiffness on the descending limb of the active-force-length curve. In addition, this experiment also highlights the slow recovery of the muscle’s force after stretching has ceased, and the phenomena of passive force enhancement after stimulation is removed. Here we will examine the 9 mm/s ramp experiment in detail because the simulations of the 3 mm/s and 27 mm/s ramp experiments produces similar stereotypical patterns (see Appendix 7 for details).

When the active lengthening experiment of Herzog and Leonard, 2002 is simulated (Figure 9A), both models produce a force transient initially (Figure 9B), but for different reasons. The VEXAT model’s transient is created when the lumped crossbridge spring (the ${\displaystyle {\tilde{k}}_{\text{o}}^{\text{X}}{\tilde{\mathcal{l}}}^{\text{X}}}$ term in Equation 15) is stretched. In contrast, the Hill model’s transient is produced, not by spring forces, but by damping produced by the force-velocity curve as shown in Equation 26.

Figure 9

Download asset Open asset

After the initial force transient, the response of the two models diverges (Figure 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 (Figure 1C). As a result, the ${\displaystyle {\mathcal{l}}^{\text{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 (Figure 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 (Figure 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 (Figure 9D). Equation 27 and Figure 2B shows the reason that damping drops during lengthening: ${\displaystyle d{\mathbf{f}}^{\text{V}}/d{v}^{\text{M}}}$, the slope of the line in Fig. Figure 2B, is quite large when the muscle is isometric and becomes quite small as the rate of lengthening increases.

After the ramp stretch is completed (at time 3.4 s in Figure 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 ${\displaystyle {\beta }_{\text{A}}^{\text{PEVK}}}$ to ${\displaystyle 71.9f_{\text{o}}^{\text{M}}/({\mathcal{l}}_{\text{o}}^{\text{M}}/s)}$ for the elastic-tendon model, and ${\displaystyle 77.7f_{\text{o}}^{\text{M}}/({\mathcal{l}}_{\text{o}}^{\text{M}}/s)}$ for the rigid-tendon model, to match the rate of force decay of the cat soleus in the data of Herzog and Leonard, 2002. The Hill model, in contrast, recovers to its isometric value quite rapidly. Since the Hill model’s force enhancement during lengthening is a function of the rate of lengthening, when the lengthening ceases, so too does the force enhancement.

Once activation is allowed to return to zero, the data of Herzog and Leonard, 2002 shows that the cat soleus continues to develop a tension that is ${\displaystyle \mathrm{\Delta }{f}_B}$ above passive levels (Figure 9B for ${\displaystyle t>8.5s}$). The force ${\displaystyle \mathrm{\Delta }{f}_B}$ is known as passive force enhancement, and is suspected to be caused by titin binding to actin (Herzog, 2019). 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 ${\displaystyle \mathrm{\Delta }{f}_B}$ below the experimental data of Herzog and Leonard (Figure 9B) after lengthening and activation have ceased.

Active lengthening beyond actin-myosin overlap

One of the great challenges that remains is to decompose how much tension is developed by titin (Figure 1C) separately from myosin (Figure 1B) in an active sarcomere. The active-lengthening experiment of Leonard et al., 2010 provides some insight into this force distribution problem because they recorded active forces both within and far beyond actin-myosin overlap. The data of Leonard et al., 2010 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 ${\displaystyle 3.38{\mathcal{l}}_{\text{o}}^{\text{M}}}$ and force of ${\displaystyle 5.14f_{\text{o}}^{\text{M}}}$, on average. In contrast, during passive lengthening myofibrils failed at a much shorter length of ${\displaystyle 2.86{\mathcal{l}}_{\text{o}}^{\text{M}}}$ with a dramatically lower tension of of ${\displaystyle 1.31f_{\text{o}}^{\text{M}}}$. To show that the extraordinary forces beyond actin-myosin overlap can be ascribed to titin, Leonard et al., 2010 repeated the experiment but deleted titin using trypsin: the titin-deleted myofibrils failed at short lengths and insignificant stresses. Using the titin model of Equation 20 (Figure 1A) as an interpretive lens, the huge forces developed during active lengthening would be created when titin is bound to actin leaving the distal segment of titin to take up all of the strain. Conversely, our titin model would produce lower forces during passive lengthening because the proximal Ig, PEVK, and distal Ig regions would all be lengthening together (Figure 3A).

Since the experiment of Leonard et al., 2010 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 8 for parameter values). To simulate a rabbit myofibril we created a force-length curve (Rassier et al., 1999) consistent with the filament lengths of rabbit skeletal muscle (Higuchi et al., 1995; 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., 2005 of rabbit psoas titin consisting of a 70–30% mix of a 3300kD and a 3400kD titin isoform (see Appendix 2.4 for fitting details and Appendix 8 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 the ‘Model’ section, because this experiment includes extreme lengths, we consider two different force-length relations for each segment of titin (Figure 10A): a linear extrapolation, and an extension that follows the WLC model. While both versions of the titin model are identical up to ${\displaystyle {\tilde{\mathcal{l}}}_{\text{toe}}^{\text{PE}}}$, beyond ${\displaystyle {\tilde{\mathcal{l}}}_{\text{toe}}^{\text{PE}}}$ the WLC model continues to develop increasingly large forces until all of the Ig domains and PEVK residues have been unfolded and the segments of titin reach a physical singularity: at this point the Ig domains and PEVK residues cannot be elongated any further without breaking molecular bonds (see Appendix 2.4 for details). Our preliminary simulations indicated that the linear titin model’s titin-actin bond was not strong enough to support large tensions, and so we increased the value of ${\displaystyle {\beta }_{\text{A}}^{\text{PEVK}}}$ from 71.9 to 975 (compare Table 1H and Appendix 8—table 1G).

Figure 10

Download asset Open asset

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 the data of Leonard et al., 2010. No adjustments were made to the active elements of the Hill model.

When the slow active stretch (0.1 µm/sarcomere/s) of the experiment of Leonard et al., 2010 is simulated, only the VEXAT model with the linear titin element can match the experimental data of Leonard et al., 2010 (Figure 10B). The Hill model cannot produce active force for lengths greater than ${\displaystyle 1.62{\mathcal{l}}_{\text{o}}^{\text{M}}}$ since the active force-length curve goes to zero (Figure 2A and Figure 10B) and the model lacks any element capable of producing force beyond this length. In contrast, the linear titin model continues to develop active force until a length of 3.38 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ is reached, at which point the titin-actin bond is pulled off the end of the actin filament and the active force is reduced to its passive value.

The WLC titin model is not able to reach the extreme lengths observed by Leonard et al., 2010. 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 1.99 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ (Figure 10B). After 1.99 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ (Figure 10B), the tension of the WLC titin model drops to its passive value but continues to increase until the contour lengths of all of the segments of titin are reached at 2.32 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. Comparing the response of the linear model to the WLC titin model two things are clear: the linear titin model more faithfully follows the data of Leonard et al., 2010, but does so with titin segment lengths that exceed the maximum contour length expected for the isoform of titin in a rabbit myofibril.

This simulation has also uncovered a surprising fact: the myofibrils in the experiments of Leonard et al., 2010 do not fail at 2.32 ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$, as would be expected by the WLC model of titin, but instead reach much greater lengths (Figure 10B). Physically, it may be possible for a rabbit myofibril to reach these lengths (without exceeding the contour lengths of the proximal Ig, PEVK, and distal Ig segments) if the bond between the distal segment of titin and myosin breaks down. This would allow the large Ig segment, that is normally bound to myosin, to uncoil and continue to develop the forces observed by Leonard et al., 2010. Unfortunately the mechanism which allowed the samples in Leonard et al.’s experiments to develop tension beyond titin’s contour length remains unknown.

Force-length and force-velocity

Although the active portion of the Hill model is embedded in Equation 16, it is not clear if the VEXAT model can still replicate force-velocity experiments of Hill, 1938 and the force-length experiments of Gordon et al., 1966. Here, we simulate both of these experiments using the cat soleus model that we have used for the simulations described in the ‘Model’ section (‘Stochastic length perturbation experiments’) and compare the results to the force-length and force-velocity curves that are used in the Hill model and in Equation 16 of the VEXAT model.

The force-velocity experiment of Hill, 1938 is simulated by activating the model, and then by changing its length to follow a shortening ramp and a lengthening ramp. During shortening experiments, the CE shortens from ${\displaystyle 1.1{\mathcal{l}}_{\text{o}}^{\text{M}}}$ to ${\displaystyle 0.9{\mathcal{l}}_{\text{o}}^{\text{M}}}$ with the measurement of active muscle force is made at ${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$. Lengthening experiments are similarly made by measuring muscle force mid-way through a ramp stretch that begins at ${\displaystyle 0.9{\mathcal{l}}_{\text{o}}^{\text{M}}}$ and ends at ${\displaystyle 1.1{\mathcal{l}}_{\text{o}}^{\text{M}}}$. When an elastic-tendon model is used, we carefully evaluate initial and terminal path lengths to accommodate for the stretch of the tendon so that the CE still shortens from ${\displaystyle 1.1{\mathcal{l}}_{\text{o}}^{\text{M}}}$ to ${\displaystyle 0.9{\mathcal{l}}_{\text{o}}^{\text{M}}}$ and lengthens from ${\displaystyle 0.9{\mathcal{l}}_{\text{o}}^{\text{M}}}$ to ${\displaystyle 1.1{\mathcal{l}}_{\text{o}}^{\text{M}}}$.

The VEXAT model produces forces that differ slightly from the ${\displaystyle {\mathbf{f}}^{\text{V}}}$ that is embedded in Equation 16 while the Hill model reproduces the curve (Figure 11). The maximum shortening 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 ${\displaystyle {\mathbf{f}}^{\text{V}}}$ curve more rapidly by decreasing ${\displaystyle {\tau }^{\text{S}}}$ this has the undesirable consequence of degrading the low-frequency response of the model when simulating the experiments of Kirsch et al., 1994 (particularly Figure 6C., and F.). These small differences can be effectively removed by scaling ${\displaystyle v_{max}^{\text{M}}}$ by ${\displaystyle {\text{s}}^{\text{V}}}$ (Figure 11A has ${\displaystyle {\text{s}}^{\text{V}}=0.95}$) to accommodate for the small decrease in force caused by the viscoelastic XE element.

Figure 11

Download asset Open asset

The force-length experiments of Gordon et al., 1966 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 the experiment of Gordon et al., 1966, we subtracted off the passive force from the active force before producing the active-force-length profile.

The simulation of Gordon et al., 1966 shows that the VEXAT model (Figure 12A, blue dots) produces a force-length profile that is shifted to the right of the Hill model (Figure 12A, red line) due to the series elasticity introduced by the XE. We can solve for the size of this rightwards shift by noting that Equation 16 will drive the ${\displaystyle {\tilde{\mathcal{l}}}^{\text{S}}}$ to a length such that the isometric force developed by the XE is equal to that of the embedded Hill model

(30) $a{\mathbf{f}}^{\text{L}}({\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}}){\tilde{k}}_{\text{o}}^{\text{X}}{\tilde{\mathcal{l}}}^{\text{X}}=a{\mathbf{f}}^{\text{L}}({\tilde{\mathcal{l}}}^{\text{S}}+{\tilde{\text{L}}}^{\text{M}})$

allowing us to solve for

(31) ${\displaystyle {\tilde{\mathcal{l}}}^{\text{X}}=\frac{1}{{\tilde{k}}_{\text{o}}^{\text{X}}}}$

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 ${\displaystyle {\mathbf{f}}^{\text{L}}}$ curve by a constant ${\displaystyle \frac{2}{{\tilde{k}}_{\text{o}}^{\text{X}}}}$. This shift, ${\displaystyle {\mathrm{\Delta }}^{\text{L}}}$, can be calibrated out of the model (Figure 12A, magenta squares) by adjusting the ${\displaystyle {\mathbf{f}}^{\text{L}}(\cdot )}$ curve so that it is ${\displaystyle \frac{2}{{\tilde{k}}_{\text{o}}^{\text{X}}}}$ to the left of its normal position. Note that all simulations described in the previous sections made use of the VEXAT model with the calibrated force-length relation and the calibrated force-velocity relation.

Figure 12

Download asset Open asset

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 (Kirsch et al., 1994; Herzog and Leonard, 2002; Leonard et al., 2010) while also retaining the ability to reproduce the force-velocity experiment of Hill, 1938 and the force-length experiment of Gordon et al., 1966. 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 spring-damper, 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 activation-dependent stiffness of muscle originates primarily from the stiffness of the attached crossbridges (Veigel et al., 1998), the origins of the activation-dependent damping observed by Kirsch et al., 1994 have not yet been established. We assumed that, since the damping observed by Kirsch et al., 1994 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 (Millard et al., 2013) given enough time (Equation 16). A recent experiment (see Figure 7 of Tomalka et al., 2021) suggests the force developed by the XE might decrease during lengthening rather than increasing as is typical of a Hill model (Millard et al., 2013). If the observations of Tomalka et al., 2021 can be replicated, the VEXAT model will need to be adjusted so that 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 Equation 16. This assumption is a departure from previous models: Hill-type models (Zajac, 1989; Millard et al., 2013) assume that the tension applied to the muscle instantaneously affects the actin-myosin sliding velocity; Huxley models (Huxley, 1957) 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., 2009 and Schappacher-Tilp et al., 2015: 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: the winding filament theory of Nishikawa et al., 2012 and the titin entanglement hypothesis of DuVall et al., 2016. 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., 1994: 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 1). 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. Although 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 Huxley and Simmons, 1971; Howard, 1997), the accuracy of the model will degrade as the number of attached crossbridges 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, Equation 16 and Equation 21, that assume myosin-actin and titin-actin interactions are temperature-invariant and scale linearly with size (${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ and ${\displaystyle f_{\text{o}}^{\text{M}}}$). In contrast, myosin-actin interactions and some titin-actin interactions are temperature-sensitive (Minajeva et al., 2002; Roots et al., 2012) and may not scale linearly with size. In the ‘Active lengthening beyond actin-myosin overlap’ section we had to adjust the active titin damping parameter, ${\displaystyle {\beta }_{\text{A}}^{\text{PEVK}}}$, to simulate the myofibril experiments of Leonard et al., 2010, perhaps because the assumptions of temperature-invariance and size-linearity were not met: the initial value for ${\displaystyle {\beta }_{\text{A}}^{\text{PEVK}}}$ came from fitting to in-situ experimental data (Herzog and Leonard, 2002) from whole muscle that was warmer (35 − 36.5°C vs 20 − 21°C) and larger (${\displaystyle {\mathcal{l}}_{\text{o}}^{\text{M}}}$ of 42.9mm vs. 10 − 15μm) than the myofibrils (Leonard et al., 2010). While the cat soleus XE and titin model parameters (Table 1G, H1) 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 (Stephenson and Wendt, 1984), and the scaling of the maximum shortening velocity (Chow and Darling, 1999). 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: the experiments of Kirsch et al., 1994 which occur over small time and length scales; and the active lengthening experiments of Herzog and Leonard, 2002 and Leonard et al., 2010 which occur over physiological and supra-physiological length scales. In contrast, we have shown in the ‘Biological benchmark simulations’ section that a Hill-type model compares poorly to biological muscle when the same set of experiments are simulated. We expect that a Huxley model (Huxley, 1957) is also likely to have difficulty reproducing the experiment of Kirsch et al., 1994 because the model lacks an active damping element. Since titin was discovered (Maruyama, 1976) long after Huxley’s model was proposed (Huxley, 1957), a Huxley model will be unable to replicate any experiment that is strongly influenced by titin such as the experiment of Leonard et al., 2010.

Although there have been several more recent muscle model formulations proposed, none have the properties to simultaneously reproduce the experiments of Kirsch et al., 1994, Herzog and Leonard, 2002, Leonard et al., 2010, Hill, 1938, and Gordon et al., 1966. Linearized impedance models (Hogan, 1985; Mussa-Ivaldi et al., 1985) can reproduce the experiments of Kirsch et al., 1994, but these models lack the nonlinear components needed to reproduce the force-length experiments of Gordon et al., 1966, and the force-velocity experiments of Hill, 1938. The models of Forcinito et al., 1998, and Tahir et al., 2018 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., 1994, in contrast, indicate that the phase shift between length and force approaches ninety degrees with increasing frequencies. Although the Hill-type models of Haeufle et al., 2014 and Günther et al., 2018 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 (Millard et al., 2013), it is likely that models of Haeufle et al., 2014 and Günther et al., 2018 will not be able to match the frequency response of biological muscle. While the model of Tamura and Saito, 2002 is one of the few models that can develop force-enhancement and force-depression (Tamura et al., 2005), it is unlikely that this model will be able to reproduce the frequency response of biological muscle because it uses spring-damping elements in series: Kirsch et al., 1994 showed that the frequency response of spring-damping elements in series poorly fits the response of biological muscle. The models of De Groote et al., 2017 and De Groote et al., 2018 introduced a short-range-stiffness element in parallel to a Hill model to capture the stiffness of biological muscle. While the formulations presented in De Groote et al., 2017 and De Groote et al., 2018 improves upon a Hill model it is unlikely to reproduce the experiement of Kirsch et al., 1994 because we have shown in the ‘Active lengthening beyond actin-myosin overlap’ section that a Hill model has a frequency response that differs from biological muscle. The muscle model of Rode et al., 2009 uses a Hill model for the CE and so we expect that this model will have the same difficulties reproducing the experiment of Kirsch et al., 1994. The model of Schappacher-Tilp et al., 2015 extends a Huxley model (Huxley, 1957) by adding a detailed titin element. Similar to a Huxley model, the model of Schappacher-Tilp et al., 2015 will likely have difficulty reproducing the experiment of Kirsch et al., 1994 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 does muscle stiffness and damping change across the force-length curve? Does stiffness and damping change with velocity? What are the physical origins of the active damping observed by Kirsch et al., 1994? 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 online (see the elife 2023 branch of https://github.com/mjhmilla/Millard2021ImpedanceMuscle; copy archived at Millard, 2024).