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. Kirsch et al. [5] were able to show that muscle behaves like a spring in parallel with a damper for small 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 [6, 7], 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 [6, 7]. 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 short-range [5] and long-range perturbations [6, 7] while also retaining the capability to reproduce the experiments of Hill [8] and Gordon et al. [9]. 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 model¹ 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 [8], mechanistic Huxley [10, 11] models in which individual elastic crossbridges are incorporated, and linearized muscle models [12, 13] which are accurate for small changes in muscle length. Kirsch et al. [5] demonstrated that, in the short-range, 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 [14, 15] 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 [8], nor Gordon et al.’s [9] 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. [16] modelled the velocity dependence of muscle using a rheological element² and an elastic rack rather than embedding the force-velocity relationship in equations directly, as is done in a typical Hill model [14, 15]. This modification allows Forcinito et al.’s [16] model to more faithfully replicate the force development of active muscle, as compared to a Hill-type model, during ramp length changes of ≈ 10%³ of the optimal CE length, and across velocities of 4 − 11% of the maximum contraction velocity⁴. Haeufle et al. [18] made use of a serial-parallel network of spring-dampers to allow their model to reproduce Hill’s force-velocity relationship [8] 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 [19] than is possible with a typical Hill model. Gü nther et al. [20] 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 [21, 22] has been more recently investigated to explain how lengthened muscle can develop active force when lengthened both within, and beyond, actin-myosin overlap [7]. 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 [23]. 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 [24] since the distal immunoglobulin (IgD) segments bind strongly to the thick filament [25]. Titin has proven to be a complex filament, varying in composition and geometry between different muscle types [26, 27], widely between species [28], and can apply activation dependent forces to actin [29]. 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 [30, 31], between titin’s PEVK region and actin [32, 33], between titin’s N2A region and actin [34], and between the PEVK-IgD regions of titin and myosin [35]. 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 [36, 37, 38, 39, 40, 41], and more recently with myosin [42].

Although the addition of a titin element to a model will result in more accurate force production during large active length changes, a titin element alone does little to affect the stiffness and damping of muscle in the short range because of its relatively low stiffness. A single cross bridge has a stiffness of 0.69 ± 0.47 pN/nm [43], which is between 22-91 times greater than the average stiffness of a single titin molecule [44] of 0.015 pN/nm from a human soleus muscle⁵. Taking into account the variability in observed stiffness of crossbridges [43], actin [45], and myosin filaments, the actin-myosin load path between the Z and M-line will be 54-233 times stiffer than the comparable titin load path assuming that only 25% of the crossbridges are attached [46]. If the middle of titin’s PEVK region rigidly attaches to actin then the actin-myosin load path is only 15-68 times stiffer than the titin load path (see Appendix A for further details). Thus it is likely that attached cross-bridges dominate active stiffness of whole muscle in the short-range [5]. Since titin-focused models have not made any changes to the modelled myosin-actin interaction beyond a Hill [14, 15] or Huxley [10, 11] model, it is unlikely that these models would be able to replicate Kirsch et al.’s experiments [5].

Although most motor control simulations [47, 48, 49, 2, 50] make use of the canonical linearized muscle model, phenomenological muscle models have also been used and modified to include stiffness. Sartori et al. [51] 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 [6]. 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 [52], 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. [53, 54] 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. [53, 54] 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. The response of the model to all perturbations within actin-myosin overlap is dominated by a viscoelastic crossbridge element that has different dynamics over short and long time-scales: over short time-scales the viscoelasticity of the lumped crossbridge dominates the response of the muscle [5], while over longer time-scales the response of the crossbridges is dominated by the force-velocity [8] and force-length [9] properties of muscle. To capture the active forces developed by muscle beyond actin-myosin overlap we added an active titin element which, similar to existing models [36, 38], features an activation-dependent⁶ 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 [55], 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. [5]. Herzog et al.’s [6] 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 [7] 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 [8] and force-length experiments of Gordon et al. [9]. Since Hill-type models are so commonly used, we also replicate all of the simulated experiments using Millard et al.’s [15] 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 [14] and makes it possible to model vastly different musculotendon units (MTUs) by simply changing the architectural and contraction properties: the maximum isometric force , the optimal CE length (at which the CE develops ), the pennation angle of the CE with respect to the tendon α_o at a CE length of , the maximum shortening velocity of the CE, and the slack length of the tendon . Many properties of sarcomeres scale with and scales with physiological cross-sectional area [56], the force-length property scales with [57], the maximum normalized shortening velocity of different CE types scales with across animals great and small [58], and that titin’s passive-force-length properties scale from single molecules to myofibrils [59, 60]

Figure 1:

The proposed model has several additional properties that we assume to scale with and inversely with : the maximum active isometric stiffness and damping , the passive forces due to the extracellular matrix (ECM), and passive forces due to titin. As crossbridge stiffness is well studied [61], we assume that muscle stiffness due to crossbridges scales such that

where is the maximum normalized stiffness. This scaling is just what would be expected when many crossbridges [61] 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. [5] is due to the intrinsic damping properties of individual crossbridges which will also scale linearly with and inversely with

where 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) cross-bridge (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 together to a single lumped crossbridge 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 [45] and myosin filaments [62] greatly exceeds that of crossbridges [43], 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 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 [14, 15], we assume that the passive elasticity of these structures scale linearly with and inversely with . Since the VEXAT model has two passive load paths (Fig. 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. [60].

As previously mentioned, there are currently several competing theories to explain how titin interacts with the other filaments in activated muscle but no definitive data to support one theory over another. While there is evidence for titin-actin interaction near titin’s N2A region [34], there is also support for a titin-actin interaction occurring near titin’s PEVK region [32, 33], and for a titin-myosin interaction near the PEVK-IgD region [35]. 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 [32, 33] in the PEVK segment; and second, there is evidence supporting an interaction at the proximal end of the PEVK segment (N2A-actin interaction) [34]. 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 Kellermayer et al.’s observations [29]: adding titin filaments and calcium slowed, but did not stop, the progression of actin filaments across an plate covered in active crossbridges (an in-vitro motility array). When activated, we assume that the amount of damping between the titin and actin scales linearly with and inversely 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 rigid tendon model given q^R, and

to evaluate the pennation angle of an elastic-tendon model. We have added a small compressive element KE (Fig. 1A to prevent the model from reaching the numerical singularity that exists in the pennation model when α → 90°. 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 ℓ^S and the length of the myosin filament L^M. The length of ℓ², 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 (ℓ¹) and the length of the IgD segment that is bound to myosin (L^IgD). 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 short time-scales, while all of the other experiments take place at much longer time-scales. Here we mathematically decouple phenomena that affect short 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) [15].

Before proceeding, a small note on notation: throughout this work we will use a 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 𝒞² continuous⁷ Bézier splines so that the model is compatible with gradient-based optimization. Normalized quantities follow a specific convention: lengths for components within the CE are normalized by the optimal CE length , the tendon is normalized by tendon slack length, velocities by the maximum shortening velocity , forces by the maximum active isometric tension , stiffness and damping by .

To evaluate the state derivative of the model, we require equations for , v¹, and v^M if the tendon is elastic. For we use of the first order activation dynamics model described in Millard et al. [15]⁸ 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 with the intention of having the model behave like a spring-damper over short periods of time, but to converge to the tension developed by a Hill-type model

over longer periods of time, where f ^L(·) is the active-forcelength curve (Fig. 2A), f ^PE(·) is the passive-force-length curve (Fig. 2A), and f ^V(·) is the force-velocity (Fig. 2B).

Figure 2:

The normalized tension developed by the VEXAT model

differs from that of a Hill model, Eqn. 14, because it has no explicit dependency on , includes four passive terms, and a lumped viscoelastic crossbridge element. The four passive terms come from the ECM element (Fig. 3A), the PEVK-IgD element (Fig. 3A and B), the compressive term (which prevents from reaching a length of 0), and a numerical damping term (where is small). The active force developed by the XE’s lumped crossbridge is scaled by the fraction of the XE that is active and attached, ), where f ^L(·) is the active-force-length relation (Fig. 2A). We evaluate f ^L using in Eqn. 15, rather than as in Eqn. 14, since actin-myosin overlap is independent of crossbridge strain. With derived, we can proceed to model the acceleration of the CE, , 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:

We set the first term of so that, over time, the CE is driven to develop the same active tension as a Hill-type model [15] (terms highlighted in blue)

where τ^S is a time constant and is the force-velocity curve (Fig. 2B). The rate of adaptation of the model’s tension, to the embedded Hill model, is set by the time constant τ^S: as τ^S is decreased the VEXAT model converges more rapidly to a Hill-type model; as τ^S is increased the active force produced by the model will look more like a springdamper. Our preliminary simulations indicate that there is a trade-off to choosing τ^S: when τ^S is large the model will not shorten rapidly enough to replicate Hill’s experiments, while if τ^S is small the low-frequency response of the model is compromised when Kirsch et al.’s [5] experiments are simulated.

The remaining two terms, and , 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, , to prevent this second-order system from unrealistically oscillating⁹. The final term , where G is a scalar gain and a_L is a low-activation threshold (a_L is 0.05 in this work), has been included as a consequence of the generalized positions we have chosen. When the CE is nearly deactivated (as a approaches a_L) this term forces and 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 ℓ^S a state variable, rather than ℓ^X. We chose ℓ^S as a state variable, rather than ℓ^X, so that the states are more equally scaled for numerical integration.

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 and and the ECM . We model titin’s elasticity as being due to two serially connected elastic segments: the first elastic segment is formed by lumping together the IgP segment and a fraction Q of the PEVK segment, while the second elastic segment is formed by lumping together the remaining (1 − Q) of the PEVK segment with the free IgD section. Our preliminary simulations of Herzog and Leonard’s active lengthening experiment [6] indicate that a Q value of a half, 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 , and 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 (Fig. 2A and Fig. 3A) [60]; that the proportion of the passive force developed by titin and the ECM is within experimental observations [60] (Fig. 3A); and that the IgP and PEVK sections show the same relative elongation as observed by Trombitás et al. [24] (Fig. 3C). Since both the passive–force-length relation and Trombitás et al.’s [24] data is at modest lengths, we also 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 the worm-like-chain (WLC) model [65]. For details, please see Appendix B.3.

When activated, we assume that some point of the PEVK segment bonds with actin through an activation-dependent damper, and that this bond forms within a specific range of CE lengths. Hisey et al. [66] observed that active lengthening produces no residual force enhancement (RFE) on the ascending limb. In the VEXAT model, the distal segment of titin, ℓ², can contribute to RFE when the titin-actin bond is formed and CE is lengthened beyond , the length at which passive force begins to develop. To incorporate Hisey et al.’s [66] observed RFE length dependence, we introduce a smooth step function

that transitions from zero to one as extends beyond . The sharpness of the transition of u^S between zero and one is controlled by R. 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, , exceeds the length of the actin filament, . To break the titin-actin bond at long CE lengths we introduce another smooth step function

The strength of the titin-actin bond also appears to vary nonlinearly with activation. Fukutani and Herzog [67] observed that the absolute RFE magnitude produced by actively lengthened fibers is similar between normal and reduced contractile force states. Since these experiments [67] 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 [67] 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 Ca²⁺ concentration [68], this is a mathematical convenience. The work of Leonard et al. [7] makes it clear that both Ca²⁺ and cross-bridge cycling are needed to allow titin to develop enhanced forces during active lengthening: no enhanced forces are observed in the presence of Ca²⁺ when cross-bridge cycling is chemically inhibited. Putting this all together, the active damping acting between the titin and actin filaments is given by the product of where is the maximum damping coefficient.

With a model of the titin-actin bond derived, we can model 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

Since Eqn. 20 is linear in , we can solve directly for it

We have added a small amount of numerical damping β^ϵ to the denominator of Eqn. 21 to prevent a singularity when the CE is deactivated.

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 of a rigid-tendon model can be evaluated directly given the path length, velocity, and the time derivatives of Eqns. 6 and 8. After v¹ is evaluated using Eqn. 21, the velocities of the remaining segments can be evaluated using the time derivatives of Eqns. 9-11.

Evaluating the CE rate of lengthening, v^M, for an elastic tendon muscle model is more involved. As is typical of lumped parameter muscle models [14, 69, 15], here we assume that difference in tension, , between the CE and the tendon

is negligible¹⁰. 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 impedances combine in series [71], 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 , is a linear scaling of the normalized tendon stiffness by U, a constant scaling coefficient. We have chosen this specific damping model because it fits the data of Netti et al. [72] and captures the structural coupling between tendon stiffness and damping (see Appendix B.1 and Fig. 15 for further details). Now that all of the terms in Eqn. 22 have been explicitly defined, we can use Eqn. 22 to solve for the remaining term, v^M, in the state derivative of the elastic-tendon model. Equation 22 becomes linear in v^M after substituting the force models described in Eqns. 23 and 15, and the kinematic model equations described in Eqns. 8, 9 and 11 (along with the time derivatives of Eqns. 8-11). After some simplification we arrive at

allowing us to evaluate the final state derivative in . During simulation the denominator of will always be finite since , and α < 90° due to the compressive element. The evaluation of in the VEXAT model is free of numerical singularities, giving it an advantage over a conventional Hill-type muscle model [15]. In addition, the VEXAT’s does not require iteration to numerically sovle a root, giving it an advantage over a singularity-free formulation of the Hill model [15]. As with previous models, initializing the model’s state is not trivial and required the derivation of a model-specific method (see Appendix C for details).

The VEXAT model introduces components and parameters that do not appear in a conventional Hill-type model model: a lumped viscoelastic XE with two parameters (and ), a two-segment active titin model with two curves and and one parameter , a tendon damping model with one parameter U (see Appendix B), and Eqn. 16 that relates the relative acceleration between actin and myosin to a force imbalance. Equation 16 introduces a time-constant (τ^S) and a damping-like term D. Fortunately, there is enough experimental data in the literature that default values for these new parameters can be evaluated and normalized using the architectural properties of the muscle (Appendix B). These default values have been derived from experiments on cats and rabbits: assuming that musculotendon dynamics look similar when normalized, these default values should be reasonable first approximations for a wide variety of mammals and muscles.

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 short-range (1-3.8% ) stochastic perturbation experiment of Kirsch et al. [5] 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 [1, 73]: co-contraction is high during initial learning, but decreases over time as a task becomes familiar. The active lengthening experiment of Herzog and Leonard [6] shows that modestly stretched (7-21% ) biological muscle has positive stiffness even on the descending limb of the active force-length curve . In contrast, a conventional Hill model [14, 15] 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. [7] unequivocally demonstrates that the CE continues to develop active forces during extreme lengthening which exceeds actin-myosin overlap. Active force development beyond actin-myosin overlap is made possible by titin, and its activation dependent interaction with actin [7]. The biological benchmark simulations conclude with a replication of the force-velocity experiments of Hill [8] and the forcelength experiments of Gordon et al. [9]. Unless otherwise noted, all simulations are performed using the same model 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 (Fig. 4). Kirsch et al. [5] applied stochastic length perturbations (Fig. 5A) to elicit force responses from the muscle (Fig. 5B) across a broad range of frequencies (4-90 Hz) and across a range of small length perturbations . Next, they applied system identification methods [71] to extract curves that defined how the gain and phase shift (between the output force of the muscle and the input length) varied as a function of frequency (Fig. 5C).

Figure 4:

Figure 5:

System identification methods make it possible to identify the mechanical system that best explains the output force waveform given the input length waveform. Using this technique, Kirsch et al. [5] were able to compare the frequency response of the specimen (See Appendix D for details) to several candidate models. They found that the pattern of gain and phase shift was most consistent with a parallel spring-damper. 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 between the two time-domain signals

Astonishingly, Kirsch et al. [5] demonstrated that a spring-damper of best fit has a VAF of between 88-99% 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)¹¹. 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 2^nd order Butterworth filter is applied to the waveform and finally the signal is scaled to the appropriate amplitude (Fig. 6). The muscle model is then activated until it develops a constant tension at a length of . 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 modelled muscles are used to evaluate the gain and phase responses of the models in the frequency-domain (See Appendix D for details). Note that the frequency resolution of our simulations are finer than the Kirsch et al.’s [5] data since we can use all 6.15 seconds of the simulated data (2048 samples at 333 Hz), while Kirsch et al. [5] were restricted to analyzing smaller time windows during the experiment (also 2048 samples at 333 Hz) in which the average output of the muscle could be considered constant.

Figure 6:

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. 7A to G), while at 90Hz the VEXAT model reaches a VAF of 91% (compare Fig. 7D to J) which is within the range of 88-99% reported by Kirsch et al. [5]. While the VEXAT model has a gain profile in the frequency-domain that closer to Kirch et al.’s data [5] than the Hill model (compare Fig. 7B to H), both models have a greater phase shift than Kirch et al.’s data [5] (compare Fig. 7C to I). The gain response of the VEXAT model follows Kirsch et al.’s data [5] closely, while the gain response of the Hill model has widely scattered coefficients due to the nonlinearity of the model’s response to even small perturbations (compare Fig. 7E to K). The phase response of the VEXAT model to the 90 Hz perturbation (Fig. 7F) 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. 7A. and D.) and the Hill model (listed on Fig. 7E. and H., respectively).

Figure 7:

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 improve slightly (compare to Fig. 7A-F to Fig. 18A-F), the response of the Hill model to the 15 Hz perturbation changes dramatically with the time-domain VAF rising from 51% to 81%. 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. 18H), 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 by using amounts of damping that exceeds that of biological muscle as indicated by Kirsch et al.’s analysis [5].

The gain and phase profiles of both models contain coefficients that do not follow a clear line, but instead are scattered. This scattering of coefficients occurs when the underlying plant being analyzed is nonlinear, even when subject to small perturbations. Much of the VEXAT model’s nonlinearities in this experiment come from the tendon tendon model (compare to Fig. 7A-F to Fig. 18A-F), since the response of the VEXAT model with a rigid tendon is less scattered. The Hill model’s nonlinearities originate from the underlying expressions for stiffness and damping of the Hill model. The stiffness of a Hill model’s CE

is heavily influenced by the partial derivative of which has a region of negative stiffness. Although is well approximated as being linear for small length changes, changes sign across . Similarly, the damping of a Hill model’s CE

also suffers from high degrees of nonlinearity for small perturbations about v^M = 0 since the slope of is positive and large when shortening, and positive and small when lengthening (Fig. 2B).

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. 8). The VEXAT model develops stiffness and damping similar to Kirsch et al.’s data when coupled with either a viscoelastic tendon (Fig. 8A & B) or a rigid tendon (Fig. 8C & D), and with a high VAF. In contrast, when the Hill model is coupled with an elastic tendon, its damping is notably larger than Kirsch et al.’s data as is the model’s stiffness at the higher tensions (Fig. 8B). This pattern changes when simulating a Hill model with a rigid tendon: the model’s stiffness is now half of Kirsch et al.’s data (Fig. 8C), while the model’s final damping coefficient is nearly three times the value measured by Kirsch et al. (Fig. 8D). The tendon model also affects the VAF of the Hill model to a large degree: the elastic-tendon Hill model has a low VAF 31%-44% (Fig. 8A & B) while the rigid-tendon Hill model has a much higher VAF of 79%. While the VEXAT model’s stiffness and damping increases linearly with activation, consistent with Kirsch et al.’s data [5], the stiffness and damping of the Hill model differs from Kirsch et al.’s data [5], ranging from small amounts of error (rigid-tendon model developing low tensions shown in Fig. 8C & D) to large amounts of error (elastic-tendon model under high tension shown in Fig. 8A & B).

Figure 8:

When the VAF of the VEXAT and Hill model is evaluated across a range of tensions (ranging from 2.5-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 VAF of 80-100% of the VEXAT model is closer to the range of 88-99% reported by Kirsch et al. [5] than the range of 32-91% of the Hill model (Fig. 9); and second, that there are systematic variations in VAF, stiffness, and damping across the different perturbation magnitudes and frequencies as shown in Tables 4 and 5 (Appendix E). While the VAF of the VEXAT model is similar whether a rigid or elastic tendon is used across all frequencies (Fig. 9A, B, and C), the Hill model’s VAF noticeably varies with tendon type: in the lower frequency range (Fig. 9A and B) the elastic-tendon Hill model has a far lower VAF than a rigid-tendon Hill model; while in the highest frequency range the VAF of the Hill model (Fig. 9C) is not as affected by the elasticity of the tendon. While the VEXAT model’s lowest VAF occurs in response to the low frequency perturbations (Fig. 9A), 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 high frequency perturbations (Fig. 9C) while the elastictendon Hill model has its lowest VAF in response to the low frequency perturbations (Fig. 9A). It is unclear if biological muscle displays systematic shifts in VAF since Kirsch et al. [5] did not report the VAF of each trial, instead indicating that a spring-damper of best-fit displayed a VAF of 88-99% across all trials and specimens.

Figure 9:

3.2 Active lengthening on the descending limb

We now turn our attention to the active lengthening in-situ experiments of Herzog and Leonard [6]. During these experiments a cat soleus was actively lengthened by modest amounts starting on the descending limb of the active-force-length curve in Fig. 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 k^M as shown in Eqn. 26 as ℓ^M extends beyond . Herzog and Leonard’s [6] experiment 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 F for details.

When Herzog and Leonard’s [6] active ramp-lengthening (Fig. 10A) experiment is simulated, both models produce a force transient initially (Fig. 10B), but for different reasons. The VEXAT model’s transient is created when the lumped crossbridge spring (the term in Eqn. 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 Eqn. 26.

Figure 10:

After the initial force transient, the response of the two models diverges (Fig. 10B): 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 ℓ² 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. 10B), the musculotendon unit develops slightly more tension using a rigid tendon because the strain of the musculotendon unit is soley borne by the CE.

In contrast, the Hill model develops less force during lengthening when it enters a small region of negative stiffness (Fig. 10B 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. 10D). Equation 26 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 isometric and becomes quite small as the rate of lengthening increases.

After the ramp stretch is completed (at time 3.4 seconds in Fig. 10B), 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 to for the elastic tendon model, and for the rigid tendon model, to match the rate of force decay of the cat soleus in Herzog and Leonard’s data [6]. 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, Herzog and Leonard’s data shows that the cat soleus continues to develop a tension that is Δf_B above passive levels (Fig. 10B for t > 8.5s). The force Δf_B is known as passive force enhancement, and is suspected to be caused by titin binding to actin [74]. 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 Δf_B below the experimental data of Herzog and Leonard (Fig. 10B) 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 [7] active-lengthening experiment provides some insight into this force distribution problem because they recorded active forces both within and far beyond actinmyosin overlap. Leonard et al.’s [7] 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 and force of , on average. In contrast, during passive lengthening myofibrils failed at a much shorter length of with a dramatically lower tension of of . To show that the extraordinary forces beyond actin-myosin overlap can be ascribed to titin, Leonard et al. [7] repeated the experiment but deleted titin using trypsin: the titin-deleted myofibrils failed at short lengths and insignificant stresses. Using the titin model of Eqn. 20 as an interpretive lens (Fig. 1A), 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 (Fig. 3A). 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 (Fig. 3A).

Since Leonard et al.’s experiment [7] was performed on skinned rabbit myofibrils and not on whole muscle, both the VEXAT and Hill models had to be adjusted prior to simulation. To simulate a rabbit myofibril we created a force-length curve consistent [75] with the actin filament lengths of rabbit skeletal muscle (1.12µm actin, 1.63µm myosin, and 0.07µm z-line width [45]), set the dimensions of titin to be consistent with Prado et al.’s [60] measurements (50 proximal Ig domains, 800 PEVK residues, and 22 distal Ig domains), used a rigid tendon of zero length, and set the pennation angle to zero. Since the evaluation of the model is done using normalized lengths and forces, we did not scale the lengths and forces of the model to reflect length and maximum tension of the fibrils tested by Leonard et al. [7].

As mentioned in Sec. 2, because this experiment includes extreme lengths, we consider two different force-length relations for each segment of titin (Fig. 11A): a linear extrapolation, and an extension that follows the WLC model. While both versions of the titin model are identical up to , beyond the WLC model continues to develop larger and larger forces until all of the Ig domains and PEVK residues have been unfolded (see Appendix B.3 for details) 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. Finally, our preliminary simulations indicated that the linear titin model’s titin-actin bond was not strong enough to support large tensions, and so we scaled the value of used to simulate an entire cat soleus by a factor of ten.

Figure 11:

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 [7]. 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 [7] only the VEXAT model with the linear titin element can match the experimental data of Leonard et al. [7] (Fig. 11B). The Hill model cannot produce active force for lengths greater than since the active force-length curve goes to zero (Fig. 2B) 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 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. [7]. 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 (Fig. 11B). After (Fig. 11B), 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 . 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. [7], 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 Leonard et al.’s [7] experiments do not fail at , as would be expected by the WLC model of titin, but instead reach much greater lengths (Fig. 2B). 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. [7]. Unfortunately the precise mechanism which allowed the samples in Leonard et al.’s experiments to develop tension beyond titin’s contour length remains unknown.

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 [8] and Gordon et al.’s [9] 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 [8] 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 to with the measurement of active muscle force is made at . Lengthening experiments are similarly made by measuring muscle force mid-way through a ramp stretch that begins at and ends at . 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 to and lengthens from to .

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. 12). 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 f ^V curve more rapidly by decreasing τ^S this has the undesirable consequence of degrading the low-frequency response of the model when Kirsch et al.’s experiments [5] (particularly Fig. 7C., and F.). These small differences can be effectively removed by scaling (Fig. 12 shows that scaling by 1.05 results in a force-velocity curve that closely follows the desired curve) to accommodate for the small decrease in force caused by the viscoelastic XE element.

Figure 12:

Gordon et al.’s [9] 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 short 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 [9] 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 [9] experiment shows that the VEXAT model (Fig. 13A, blue dots) produces a force-length profile that is shifted to the right of the Hill model (Fig. 13A, red line) due to the series elasticity introduced by the XE. We can solve for the size of this right-wards shift by noting that Eqn. 16 will drive the to a length such that the isometric force developed by the XE is equal to that of the embedded Hill model

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 . This shift can be calibrated out of the model (Fig. 13 upper plot, magenta squares) adjusting the f ^L(·) curve so that it is 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 13:

Figure 14:

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, 6, 7] while also retaining the ability to reproduce Hill’s force-velocity [8] experiment and Gordon et al.’s [9] 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 model that over short 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 [43], 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 [15] given enough time (Eqn. 16). A recent experiment of Tomalka et al. [76] suggests the force developed by the XE might decrease during lengthening rather than increasing as is typical of a Hill model [15]. If Tomalka et al.’s [76] 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 [14, 15] assume that the tension applied to the muscle instantaneously affects the actin-myosin sliding velocity; Huxley models [10] 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 assumes that some part of the PEVK segment interacts with actin through an activation dependent damping force: the nature of this interaction is not clear, and many different mechanisms have been proposed [36, 37, 38, 42]. The active damping element we have proposed is similar to Rode et al.’s [36] sticky-spring model, but instead acts at a single location within the PEVK segment like the model of Schappacher-Tilp et al. [38]. What distinguishes our titin model is that it requires only a single state and its curves have been carefully constructed to satisfy three constraints: that the passive force-length curve of the ECM and titin together can be fit to measurements, ratio of the ECM to titin is within the range measured by Prado et al. [60], and the relative strains of the proximal Ig and PEVK sections are consistent with Trombitás et al.’s measurements [24].

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; the active lengthening experiments of Herzog and Leonard [6] and Leonard et al. [7] which occur over physiological to 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 [10] 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 [21] long after Huxley’s model was proposed [10], a Huxley model will be unable to replicate any experiment that is strongly influenced by titin such as Leonard et al.’s experiment [7].

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 [6], Leonard et al. [7], Hill [8], and Gordon et al. [9]. Linearized impedance models [12, 13] can reproduce Kirsch et al.’s experiments [5], these models lack the nonlinear components needed to reproduce Gordon et al.’s force-length curve [9] and Hill’s force-velocity curve [8]. The models of Forcinito et al. [16], Haeufle et al. [18], Gü nther et al. [20], and Tahir et al. [39] all have a structure that places a damping element in series with a spring: Kirsch et al. [5] explicitly demonstrated that this type of model fails to reproduce the gain and phase response of biological muscle. De Groote et al. [53, 54] 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 [53, 54] 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 [36] 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 [38] extends a Huxley model [10] by adding a detailed titin element. Similar to a Huxley model, Schappacher-Tilp et al.’s model [38] 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 force-length 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 online¹².

Acknowledgements

Financial support from Deutsche Forschungsgemeinschaft Grant No. MI 2109/1-1, the Lighthouse Initiative Geriatronics by StMWi Bayern (Project X, grant no. 5140951), and NSERC of Canada is gratefully acknowledged.

A The stiffness of the actin-myosin and titin load paths

A single half-myosin can connect to the surrounding six actin filaments through nearly 100 crossbridges. A 0.800 µm half-myosin has a pair crossbridges over 0.700µm of its length every 14.3 nm which amounts to 97.9 per half-myosin [77]. Assuming a duty cycle of 25% (values between 5-90% have been reported [46]), at full actin-myosin overlap there will be 24.5 attached crossbridges with a total stiffness of 5.4 − 28.4 pN/nm (24.5 × 0.69 ± 0.47 pN/nm [43]).

At full overlap, the Z-line is 1 actin filament length L ^A (1.27 µm in human [75]) from the M-line and the average point of crossbridge attachment is in the middle of the half-myosin 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 L ^A − 0.45 µm from the Z-line. A single actin filament has a stiffness of 46-68 pN/nm [45]. Since stiffness scales inversely with length, the section between the Z-line and the average attachment point has a stiffness of 83.6-123 pN/nm. Assuming that load is evenly distributed, each actin filament contributes a third of its stiffness to each of the three half-sarcomeres that it borders. Since there are six actin filaments surrounding each myosin filament the total stiffness of the actin filaments for a half-sarcomere comes to 167-246 pN/nm. Myosin has a similar stiffness as a single actin filament [62], with the section between the average attachment point and the M-line having a stiffness of 81.8-121 pN/nm. The final active stiffness of half-sarcomere comes from the series connection of actin, 24.5 crossbridges, and myosin is 4.9-21.0 pN/nm.

The stiffness of the actin-myosin load path is between 54-233 times stiffer than the average stiffness of 0.09 nN/nm generated by the six free titin filaments from skeletal muscle (each with an average stiffness of 0.015 pN/nm [44]) that span each half-sarcomere. The stiffness of titin in skeletal muscle can change dramatically during activation. If we assume that the middle of the PEVK segment is rigidly fixed to actin during activation, the stiffness of each titin filament is 3.4 times higher since only half the PEVK region and the distal Ig segment is free to stretch. Despite this large increase in stiffness, the actin-myosin load path remains 15-68 times stiffer than the titin load path.

B Model Fitting

Many of the experiments simulated in this work [5, 6] have been performed using a cat soleus muscle. While we have been able to extract some architectural parameters directly from the experiments we simulate ( and from [6]), we have had to rely on the literature mentioned in Table 1 for the remaining parameters. The remaining properties of the model can be solved by first building a viscoelastic damping model of the tendon; next, by solving for the intrinsic stiffness and damping properties of the CE; and finally, by fitting the passive curves and to simultaneously fit the passive force-length curve recorded by Herzog and Leonard [6], using a mixture of tension from titin and the ECM that is consistent with Prado et al.’s data [60], all while maintaining the geometric relationship between f ^IgP and f ^PEVK as measured by Trombitás et al. [24].

Table 1:

B.1 Fitting the tendon’s stiffness and damping

Similar to previous work [15], we model the force-length relation of the tendon using a quintic Bézier spline (Fig. 15A) that begins at (where is tendon length normalized by , and is tension normalized by ), ends at with a normalized stiffness of , and uses the constants and (given from Scott and Loeb [79], is thus 4.58%). Using the experimental data of Netti et al. [72] we have also constructed a curve to evaluate the damping coefficient of the tendon. The normalized tendon stiffness (termed storage modulus by Netti et al. [72]) and normalized tendon damping (termed loss modulus by Netti et al. [72]) both have a similar shape as the tendon is stretched from slack to (Fig. 15B and C). The similarity in shape is likely not a coincidence.

Figure 15:

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 MTUis structured with a CE in series with a tendon, the compliance of the MTUis given by

where is the stiffness of the CE in the direction of the tendon. We can calculate k^M directly by fitting a line to the stiffness vs tension plot that appears in Figure 12 of Kirsch et al. [5] (0.8mm, 0-35 Hz perturbation) and resulting in k^M =2.47 N/mm at a nominal force of 5N. Here we use a nominal tension of 5N so that we can later compare our model to the 5N frequency response reported in Figure 3 of Kirsch et al. [5]. Since Kirsch et al. [5] did not report the architectural properties of the experimental specimens, we assume that the architectural properties of the cat used in Kirsch et al.’s experiments are similar to the properties listed in Table 1. We evaluate the stiffness of the tendon model by stretching it until it develops the nominal tension of Kirsch et al.’s Figure 3 data (5N), and then compute its derivative which amounts to k^T =16.9 N/mm. Finally, using Eqn. 31 we can solve for a value of . Since the inverse of damping adds for damping elements in series

we can use a similar procedure to evaluate , the damping of the CE along the tendon. The value of β^M that best fits the damping vs. tension plot that appears in Figure 12 of Kirsch et al. [5] at a nominal tension of 5N is 0.0198 Ns/mm. The tendon damping model we have just constructed develops 0.697 Ns/mm at a nominal load of 5N. Using Eqn. 32, we arrive at . Due to the pennation model, the stiffness and damping values of the CE differ from those along the tendon.

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 by solving Eqn. 34 for k^M and substituting the values of Eqns. 37, and 39. In this case, the values of k^M (4.46 N/mm) and (4.46 N/mm) are the same to three significant figures.

We can use a similar process to transform into using the pennation model by noting that

which expands to a much smaller expression

than Eqn. 34 since α does not depend on v^M, and thus ∂α/∂v^M = 0. By taking a time derivative of Eqn. 38 we arrive at

which allows us to solve for

By recognizing that

and using Eqns. 42 and 44 we can evaluate β^M in terms of Similar to k^M, the value of β^M (0.020 Ns/mm) is close to (0.020 Ns/mm). When this same procedure is applied to the stiffness and damping coefficients extracted from the gain and phase profiles from Figure 3 of Kirsch et al. [5], the values of k^M and β^M differ (4.46 N/mm and 0.0089 Ns/mm) from the results produced using the data of Figure 12 (2.90 N/mm and 0.020 Ns/mm). Likely these differences arise because we have been forced to use a fixed maximum isometric force for all specimens when, in reality, this property varies substantially. We now turn our attention to fitting the titin and ECM elements, since we cannot determine how much of k^M and β^M are due to the XE until the titin and ECM elements have been fitted.

B.3 Fitting the force-length curves of titin’s segments

The nonlinear force-length curves used to describe titin ( and in series), and the ECM must satisfy three conditions: the total force-length curve produced by the sum of the ECM and titin must match the observed passive-force-length relation[6]; the proportion of titin’s contribution relative to the ECM must be within measured bounds [60]; and finally the stiffness of the must be a linear scaling of to match the observations of Trombitás et al. [24].

First, we fit the passive force-length curve to the data of Herzog and Leonard [6] to serve as a reference. The curve f ^PE begins at the normalized length and force coordinates of with a slope of 0, ends at with a slope of , and is linearly extrapolated outside of this region.

We solve for the and such that

the squared differences between f ^PE and the passive force-length data of Herzog and Leonard [6] (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. [60] measured across 5 different rabbit skeletal muscles¹³. The remaining fraction, 1 − P, of the force-length curve comes from titin.

In mammalian skeletal muscle, titin has three elastic segments [60] 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. [24] 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 [24] 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 [24] 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 [34, 29] and crossbridge attachment [7].

We begin by digitizing the data of Figure 5 of Trombitás et al. [24] 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 . The coefficients that best fit the data from Trom-bitás et al. [24] appear in Table 2.

Table 2:

Table 3:

These functions can be scaled to fit a titin filament of a differing geometry. Many of the experiments simulated in this work used a cat soleus. Although the lengths of the segments in a cat soleus titin filament 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 [24]) since both muscles contain predominately slow-twich fibers: slow twitch fibers tend to express longer, more compliant titin filaments [60]. 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. 51-53 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. [7] has a known titin geometry (50 proximal Ig domains, 800 PEVK residues, and 22 distal Ig domains [60]) which differs substantially from the isoform of titin expressed in the human soleus. To create the rabbit psoas titin length functions , and we begin by scaling the human soleus PEVK length function by the relative proportion of PEVK residues of 800/2174. The length of the two Ig segments

is what remains from the half-sarcomere once the rigid lengths of titin (0.100 µm for L^T12 and 0.8150 µm for L^IgD [45]) and the PEVK segment length have been subtracted away. The function that describes and can then be formed by scaling by the proportion of Ig domains in each segment

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 [24] 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 [24]). Each domain appears to be a simple nonlinear spring until it unfolds and elongates by nearly 25 nm in the process [80]. Unfolding events appear to happen individually during lengthening experiments [80], 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 [38] in total.

Active and passive lengthening experiments at the sarcomerelevel hide the complexity that is apparent when studying individual titin filaments. The experiments of Leonard et al. [7] show that the sarcomeres in a filament (from rabbit psoas) mechanically fail when stretched passively to an average length of but can reach when actively lengthened. Leonard et al. [7] showed that titin was the filament bearing these large forces since the sarcomeres were incapable of developing active or passive tension when the experiment was repeated after the titin filaments were chemically cut. It is worth noting that the forces measured by Leonard et al. [7] contain none of the complex saw-tooth pattern indicative of unfolding events even though 72 of these events would occur as each proximal and distal Ig domain fully unfolded and reached its maximal length¹⁴. Although we cannot be sure how many unfolding events occurred during Leonard et al.’s experiments [7], due to sarcomere non-homogeneity [81], it is likely that many Ig unfolding events occurred since the average sarcomere length at failure was longer than the maximum length of that would be predicted from the geometry of rabbit psoas titin¹⁵.

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 [6, 7], but will lack the nonlinear saw-tooth force-length profile that is measured when individual titin filaments are lengthened [80]. To have the broadest possible application, we will fit titin’s force-length curves to provide reasonable results for both moderate [6] and large active stretches [7]. 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 and no unfolding events occur as was done by Trombitás et al. [65]. Unfortunately, we cannot accurately capture the force development of the Ig segments at both short and long lengths without introducing many extra states to the model so that Ig domain unfolding and folding events can be simulated.

To ensure that the serially connected force-length curves of and closely reproduce , we are going to use affine transformations of f ^PE to describe and . The total stiffness of the half-sarcomere titin model is given by

which is formed by the series connection of and

Since each of titin’s segments is exposed to the same tension in Trombitás et al.’s experiment [24] the slopes of the lines that Eqns. 51-53 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 curves are beyond the toe region the stiffness is a constant and is given by

and

Dividing Eqn. 61 by 62 eliminates the unknown and results in an expression that relates the ratio of the terminal linear stiffness of and

to the relative changes in Eqns. 59 and 60. Substituting Eqns. 63, and 58 into Eqn. 57 yields the expression

which can be simplified to

and this expression can be evaluated using the terminal stiffness of titin and the coefficients listed in Table 2. Substituting Eqn. 65 into Eqn. 63 yields

The curves and can now be formed by scaling and shifting the total force-length curve of titin (1 − P)f ^PE. By construction, titin’s force-length curve develops a tension of (1 − P), and has reached its terminal stiffness, when the CE reaches a length such that . Using Eqns. 51-53 and the appropriate coefficients in Table 2 we can evaluate the normalized length developed by the ℓ¹ segment

and ℓ² segment

at a CE length of . The curve is formed by shifting and scaling the (1 − P)f ^PE curve so that it develops a normalized tension of (1 − P) and a stiffness of at a length of . Similarly, the curve is made by shifting and scaling the (1 − P)f ^PE curve to develop a normalized tension of (1 − P) and a stiffness of at a length of .

By construction, the spring network formed by the , and curves follows the fitted f ^PE curve (Fig. 3A) such that the ECM curve makes up 56% of the contribution. When the CE is active and is effectively fixed in place, the distal segment of titin contributes higher forces since undergoes higher strains (Fig. 3A). Finally, when the experiment of Trombitás et al. [24] are simulated the movements of the IgP/PEVK and PEVK/IgD boundaries in the titin model closely follow the data (Fig. 3C).

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 [24]). Scenarios in which the CE reaches extremely long lengths, such as during injury or during Leonard et al.’s experiment [7], require fitting titin’s force-length curve beyond the typical ranges observed in-vivo. The WLC model has been used successfully to model the force-length relation of individual titin segments [65] at extreme lengths. In this work, we consider two different extensions to and : a linear extrapolation, and the WLC model. Since the fitted f ^PE curve is linearly extrapolated, so too are the , and curves by default. Applying the WLC to our titin curves requires a bit more effort.

We have modified the WLC to include a slack length so that the WLC model can made to be continuous with the various force-length relations we have already derived. The normalized force developed by our WLC model is given by

where B is a scaling factor and the normalized segment length is defined as

where is the slack length, and is the contour length of the segment. To extend the curve to follow the WLC model, we first note the normalized contour length of the ℓ¹ segment

by counting the number of proximal Ig domains (N^IgP), the number of PEVK residues (QN^PEVK) associated with ℓ¹ and by scaling each by the maximum contour length of each Ig domain (25nm [80]), and each PEVK residue (between 0.32 [65] and 0.38 nm [82] 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 is given by

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. 69 so that and are continuous with each respective WLC extrapolation. However, we do not use the WLC model directly because it contains a numerical singularity which is problematic during numerical simulation. Instead, we add an additional Bézier segment to fit the WLC extension that spans between forces of (1 − P) and twice the normalized failure force noted by Leonard et al. [7]. To fit the shape of the final Bézier segment, we adjust the locations of the internal control points to minimize the squared differences between the modified WLC model and the final Bézier curve (Fig. 11A). The final result is a set of curves , and which, between forces 0 and (1 − P), will reproduce f ^PE, Trombitás et al.’s measurements [24], and do so with a reasonable titin-ECM balance [60]. For forces beyond (1 − P), the curve will follow the segment-specific WLC model up to twice the expected failure tension noted by Leonard et al. [7].

B.4 Fitting the XE’s Impedance

With the passive curves established, we can return to the problem of identifying the normalized maximum stiffness and damping of the lumped XE element. Just prior to discussing titin, we had evaluated the impedance of the cat soleus CE in Kirsch et al.’s [5] Figure 12 to be k^M =2.90 N/mm and β^M =0.020 Ns/mm at a nominal active tension of 5N. The normalized stiffness k^M can be found by taking the partial derivative of Eqn. 15 with respect to

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. 75 to

and solve for

When using to the data from Figure 12 in Kirch et al. [5], we end up with for the elastic tendon model, and for the rigid tendon model. When this procedure is repeated for Figure 3 of Kirsch et al. [5] (from a different specimen) we are left with for the elastic tendon model and for the rigid tendon model. The value for is much larger than k^M, in this case, because the a needed to generate 5N is only 0.229. Similarly, we can form the expression for the normalized damping of the CE by taking the partial derivative of Eqn. 15 with respect to

As with k^M, the expression for can be reduced to

which evaluates to for both the elastic and rigid tendon models using Kirsch et al.’s [5] Figure 12 data. The damping coefficients of the elastic and rigid tendon models is similar because the damping coefficient of the musculotendon is dominated by the damping coefficient of CE. When the data from Kirsch et al.’s [5] Figure 3 is used, the damping coefficients of the elastic and rigid tendon models are and respectively.

The dimensionless parameters , and can be used to approximate the properties of other MTUs given and . The stiffness and damping of the lumped crossbridge element will scale linearly with and inversely with provided the impedance properties of individual crossbridges, and the maximum number of crossbridges attached per unit length is similar between a feline’s skeletal muscle sarcomeres and those of the target MTU. The dimensionless titin parameters can be scaled to fit other MTUs provided these values are similar: P, A ^IgP, A ^PEVK, A ^IgD, and Q. The assumption that these parameters remain similar across different MTUs, even within the same animal [60], is crude at best. Prado et al. [60] showed that P can vary between 43%-76% across different skeletal muscles of the same rabbit. In addition, Prado et al. [60] also found that the overall stiffness of titin can vary by a factor of approximately 3.4 (see Fig. 5B [60]) and weakly correlates with the fiber type of the myofibril: fast twitch fibers tend to express a shorter (stiffer) titin, while slow twitch fibers tend to express a longer more compliant titin. As data on P is challenging to obtain, and the variation of P is large, we have chosen to use the average value of P reported by Prado et al. [60]. Since most of our simulated tests are applied to slow-twitch muscles, we have fit the stiffness of titin’s segments to the data of Trombitás et al. [24] which comes from human soleus muscle, which is often dominated by slow twitch fibers.

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 ℓ¹ 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 and using a nested bisection search that looks for values which approximately satisfies Eqn. 22, have low values of from Eqn. 16, and begins with balanced forces between the two segment titin model in Eqn. 20.

In the outer loop we iterate over values of . Given a, ℓ^P, v^P, and a candidate value of , we can immediately solve for α and ℓ^T using the pennation model. We can numerically solve for the value of another state, ℓ¹, using the kinematic relationship between ℓ^M and ℓ¹ and by assuming that the two titin segments are in a force equilibrium

Next, we iterate over values of between 0 and v^P cos α (we ignore solutions in which the sign of v^M and v^T differ) to find the value of that best satisfies Eqn. 22. Prior to evaluating Eqn. 22, we need to set both and . Here we choose a value for that will ensure that the XE is not producing transient forces

and we use fixed-point iteration to solve for such that Eqn. 16 evaluates to zero. Now the value of can be directly evaluated using the candidate value of , the first derivative of Eqn. 9, and the fact that we have set to zero. Finally, the error of this specific combination of and is evaluated using Eqn. 22, where the best solution leads to the lowest absolute value for of in Eqn. 22. If a rigid tendon model is being initialized the procedure is simpler because the inner loop iterating over is unnecessary: given v^P and are zero, the velocities and can be directly solved using the first derivative of Eqn. 9. While in principle any root solving method can be used to solve this problem, we have chosen to use the bisection method to avoid local minima.

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 and nominal force . If we approximate the muscle’s response as a linear time invariant transformation h(t) we can express

where ⋆ is the convolution operator. Each of these signals can be transformed into the frequency-domain [83] by taking the Fourier transform ℱ (·) of the time-domain signal, which produces a complex signal. Since convolution in the time-domain corresponds to multiplication in the frequency-domain, we have

In Eqn. 85 we are interested in solving for H(s). While it might be tempting to evaluate H(s) as

this is likely to result in numerical singularities, because Y (s) is only approximated by H(s)X(s): although X(s) goes to zero at high frequencies by construction (because it is low-pass filtered as described in Sec. 3.1), Y (s) may not. To evaluate H(s) in a numerically stable manner, we can multiply both sides of Eqn. 85 by X(s)

to zero the higher frequency portions of Y (s) that are nonzero as X(s) tends to zero. The signal G_xy = (X(s)Y (s)) = (Y (s)X(s)) and G_xx = (X(s)X(s)) can be formed by noting that

and

allowing us to evaluate

The gain of H(s) is given the magnitude of H(s),

while the phase is given by

where ℝ(H(s)) and 𝕀(H(s)) are the real and imaginary parts of H(s) respectively.

E Simulation summary data of Kirsch et al

Table 4:

Table 5:

F Simulations of active lengthening

Figure 16:

Figure 17:

G Supplementary Plots

Figure 18: