A three filament mechanistic model of musculotendon force and impedance

The force developed by actively lengthened muscle depends on different structures across different scales of lengthening. For small perturbations, the active response of muscle is well captured by a linear-time-invariant (LTI) system: a stiff spring in parallel with a light damper. The force response of muscle to longer stretches is better represented by a compliant spring that can fix its end when activated. Experimental work has shown that the stiffness and damping (impedance) of muscle in response to small perturbations is of fundamental importance to motor learning and mechanical stability, while the huge forces developed during long active stretches are critical for simulating and predicting injury. Outside of motor learning and injury, muscle is actively lengthened as a part of nearly all terrestrial locomotion. Despite the functional importance of impedance and active lengthening, no single muscle model has all these mechanical properties. In this work, we present the viscoelastic-crossbridge active-titin (VEXAT) model that can replicate the response of muscle to length changes great and small. To evaluate the VEXAT model, we compare its response to biological muscle by simulating experiments that measure the impedance of muscle, and the forces developed during long active stretches. In addition, we have also compared the responses of the VEXAT model to a popular Hill-type muscle model. The VEXAT model more accurately captures the impedance of biological muscle and its responses to long active stretches than a Hill-type model and can still reproduce the force-velocity and force-length relations of muscle. While the comparison between the VEXAT model and biological muscle is favorable, there are some phenomena that can be improved: the low frequency phase response of the model, and a mechanism to support passive force enhancement.

1 Introduction The stiffness and damping of muscle are properties of fun-2 damental importance for motor control, and the accurate 3 simulation of muscle force.The central nervous system 4 (CNS) exploits the activation-dependent stiffness and damp-5 ing (impedance) of muscle when learning new movements 6 [1], and when moving in unstable [2] or noisy environments 7 [3].Reaching experiments using haptic manipulanda show 8 that the CNS uses co-contraction to increase the stiffness 9 of the arm when perturbed by an unstable force field [4].10 With time and repetition, the force field becomes learned 11 and co-contraction is reduced [1].

12
The force response of muscle is not uniform, but varies 13 with both the length and time of perturbation.Under con-14 stant activation and at a consistent nominal length, Kirsch et 15 al. [5] were able to show that muscle behaves like a linear-16 time-invariant (LTI) system in response to small 1 perturba-17 tions: a spring-damper of best fit captured over 90% of the 18 observed variation in muscle force for small perturbations 19 (1-3.8%optimal length) over a wide range of bandwidths 20 .When active muscle is stretched appreciably, titin 21 can develop enormous forces [7], [8], which may prevent 22 further lengthening and injury.The stiffness that best cap-23 tures the response of muscle to the small perturbations of 24 Kirsch et al. [5] is far greater than the stiffness that best 25 captures the response of muscle to large perturbations [7], 26 [8].Since everyday movements are often accompanied by 27 both large and small kinematic perturbations, it is important 28 to accurately capture these two processes.their model to reproduce Hill's force-velocity relationship [9] mechanistically rather than embedding the experimental curve directly in their model.This modification allowed Haeufle et al.'s model to simulate high speed reaching movements that agree more closely with experimental data [21] than is possible with a typical Hill model.Günther et al. [22] evaluated how accurately a variety of spring-damper models were able to reproduce the microscopic increases in crossbridge force in response to small length changes.While each of these models improves upon the force response of the Hill model to ramp length changes, none are likely to reproduce Kirsch et al.'s experiment [5] because the linearized versions of these models lead to a serial, rather than a parallel, connection of a spring and a damper: Kirsch et al. [5] specifically showed (see Figure 3 of [5]) that a serial connection of a spring-damper fails to reproduce the phase shift between force and length present in their experimental data.Titin [23], [24] has been more recently investigated to explain how lengthened muscle can develop active force when lengthened both within, and beyond, actin-myosin overlap [8].Titin is a gigantic multi-segmented protein that spans a half-sarcomere, attaching to the Z-line at one end and the middle of the thick filament at the other end [25].In skeletal muscle, the two sections nearest to the Z-line, the proximal immunoglobulin (IgP) segment and the PEVK segment -rich in the amino acids proline (P), glutamate (E), valine (V) and lysine (K) -are the most compliant [26] since the distal immunoglobulin (IgD) segments bind strongly to the thick filament [27].Titin has proven to be a complex filament, varying in composition and geometry between different muscle types [28], [29], widely between species [30], and can apply activation dependent forces to actin [31].It has proven challenging to determine which interactions dominate between the various segments of titin and the other filaments in a sarcomere.Experimental observations have reported titin-actin interactions at myosin-actin binding sites [32], [33], between titin's PEVK region and actin [34], [35], between titin's N2A region and actin [36], and between the PEVK-IgD regions of titin and myosin [37].This large variety of experimental observations has led to a correspondingly large number of proposed hypotheses and models, most of which involve titin interacting with actin [38]- [43], and more recently with myosin [44].
The addition of a titin element to a model will result in more accurate force production during large active length changes, but does not affect the stiffness and damping of muscle at modest sarcomere lengths because of titin's relatively low stiffness.At sarcomere lengths of Since titin-focused models have not made any changes to the modeled myosin-actin interaction beyond a Hill [16], [17] or Huxley [11], [12] model, it is unlikely that these models would be able to replicate Kirsch et al.'s experiments [5].
Although most motor control simulations [2], [45]- [48] make use of the canonical linearized muscle model, phenomenological muscle models have also been used and modified to include stiffness.Sartori et al. [49] modeled muscle stiffness by evaluating the partial derivative of the force developed by a Hill-type muscle model with respect to the contractile element (CE) length.Although this approach is mathematically correct, the resulting stiffness is heavily influenced by the shape of the force-length curve and can lead to inaccurate results: at the optimal CE length this approach would predict an active muscle stiffness of zero since the slope of the force-length curve is zero; on the descending limb this approach would predict a negative active muscle stiffness since the slope of the force-length curve is negative.In contrast, CE stiffness is large and positive near the optimal length [5], and there is no evidence for negative stiffness on the descending limb of the force-length curve [7].Although the stiffness of the CE can be kept positive by shifting the passive force-length curve, which is at times used in finite-element-models of muscle [50], this introduces a new problem: the resulting passive CE stiffness cannot be lowered to match a more flexible muscle.In contrast, De Groote et al. [51], [52] modeled short-range-stiffness using a stiff spring in parallel with the active force element of a Hill-type muscle model.While the approach of de Groote et al. [51], [52] likely does improve the response of a Hill-type muscle model for small perturbations, there are several drawbacks: the short-range-stiffness of the muscle sharply goes to zero outside of the specified range whereas in reality the stiffness is only reduced [5] (see Fig. 9A); the damping of the canonical Hill-model has been left unchanged and likely differs substantially from biological muscle [5].
In this work, we propose a model that can capture the force development of muscle to perturbations that vary in size and timescale, and yet is described using only a few states making it well suited for large-scale simulations. 1 When active, the response of the model to perturbations 2 within actin-myosin overlap is dominated by a viscoelas-3 tic crossbridge element that has different dynamics across 4 time-scales: over brief time-scales the viscoelasticity of the 5 lumped crossbridge dominates the response of the muscle 6 [5], while over longer time-scales the force-velocity [9] 7 and force-length [10] properties of muscle dominate.To 8 capture the active forces developed by muscle beyond actin-9 myosin overlap we added an active titin element which, 10 similar to existing models [38], [40], features an activation-11 dependent 6 interaction between titin and actin.To ensure 12 that the various parts of the model are bounded by reality, 13 we have estimated the physical properties of the viscoelastic 14 crossbridge element as well as the active titin element using 15 data from the literature. 16 While our main focus is to develop a more accurate mus-17 cle model, we would like the model to be well suited to 18 simulating systems that contain tens to hundreds of mus-19 cles.Although Huxley models have been used to simulate 20 whole-body movements such as jumping [53], the memory 21 and processing requirements associated with simulating a 22 single muscle with thousands of states is high.Instead of 23 modeling the force development of individual crossbridges, 24 we lump all of the crossbridges in a muscle together so that 25 we have a small number of states to simulate per muscle.
To evaluate the proposed model, we compare simulations 27 of experiments to original data.We examine the response 28 of active muscle to small perturbations over a wide band-29 width by simulating the stochastic perturbation experiments 30 of Kirsch et al. [5].Herzog et al.'s [7] active-lengthening 31 experiments are used to evaluate the response of the model 32 when it is actively lengthened within actin-myosin over-33 lap.Next, we use Leonard et al.'s [8] active lengthening 34 experiments to see how the model compares to reality when 35 it is actively lengthened beyond actin-myosin overlap.In 36 addition, we examine how well the model can reproduce the 37 force-velocity experiments of Hill [9] and force-length ex-38 periments of Gordon et al. [10].Since Hill-type models are 39 so commonly used, we also replicate all of the simulated 40 experiments using Millard et al.'s [17] Hill-type muscle 41 model to make the differences between these two types of 42 models clear.We begin by treating whole muscle as a scaled half-1 sarcomere that is pennated at an angle α with respect to a 2 6 Although activation normally refers to the presence of Ca 2+ ions in the sarcomere, Ca 2+ ions alone are insufficient to cause titin to develop enhanced lengthening forces.In addition, crossbridge attachment appears to be necessary: when crossbridge attachment is inhibited titin is not able to develop enhanced forces in the presence of Ca 2+ during lengthening [8].).Titin is modeled as two springs of length ℓ 1 and ℓ 2 in series with the rigid segments L T12 and L IgD .Viscous forces act between the titin and actin in proportion to the activation of the muscle (C.), which reduces to negligible values in a purely passive muscle (D.).
We modeled actin and myosin as rigid elements; the XE, titin, and the tendon as viscoelastic elements; and the ECM as an elastic element.
tendon (Fig. 1A).The assumption that mechanical prop- erties scale with size is commonly used when modeling 4 muscle [16] and makes it possible to model vastly differ- where k X o is the maximum normalized stiffness.This scal-

26
ing is just what would be expected when many crossbridges

27
[59] act in parallel across the cross-sectional area of the 28 muscle, and act in series along the length of the muscle.

29
Although the intrinsic damping properties of crossbridges where β X o is the maximum normalized damping.For the 36 remainder of the paper, we refer to the proposed model as Similarly, we have lumped the six titin filaments per halfsarcomere (Fig. 1A) together to further reduce the number of states needed to simulate this model.
The addition of a titin filament to the model introduces an additional active load-path (Fig. 1C) and an additional passive load-path (Fig. 1D).As is typical [16], [17], we assume that the passive elasticity of these structures scale As previously mentioned, there are several theories to explain how titin interacts with the other filaments in activated muscle.While there is evidence for titin-actin interaction near titin's N2A region [36], there is also support for a titinactin interaction occurring near titin's PEVK region [34], [35], and for a titin-myosin interaction near the PEVK-IgD region [37].For the purposes of our model, we will assume a titin-actin interaction because current evidence weighs more heavily towards a titin-actin interaction than a titinmyosin interaction.Next, we assume that the titin-actin interaction takes place somewhere in the PEVK segment for two reasons: first, there is evidence for a titin-actin interaction [34], [35] in the PEVK segment; and second, there is evidence supporting an interaction at the proximal end of the PEVK segment (N2A-actin interaction) [36].We have left the point within the PEVK segment that attaches to actin as a free variable since there is some uncertainty about what part of the PEVK segment interacts with actin.
The nature of the mechanical interaction between titin and the other filaments in an active sarcomere remains uncertain.Here we assume that this interaction is not a rigid attachment, but instead is an activation dependent damping to be consistent with the observations of Kellermayer and Granzier [31] and Dutta et al. [36]: adding titin filaments and calcium slowed, but did not stop, the progression of actin filaments across a plate covered in active crossbridges (an in-vitro motility array).When activated, we assume that the amount of damping between the titin and actin scales linearly with f M o and inversely with ℓ M o .
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 Given these generalized positions, the path length ℓ P , and pennation model to evaluate the pennation angle of a CE with a rigid-tendon, and to evaluate the pennation angle of a CE with an elastic-8 tendon.We have added a small compressive element KE

9
(Fig. 1A) to prevent the model from reaching the numerical 10 singularity that exists as lM approaches lM min , the length at 11 which the pennation model when α → 90 • .The tendon 12 length of an elastic-tendon model is the difference between the 14 path length and the CE length along the tendon.The length 15 of the XE is the difference between the half-sarcomere length and the 17 sum of the average point of attachment ℓ S and the length 18 of the myosin filament L M .The length of ℓ 2 , the lumped 19 PEVK-IgD segment, is the difference between the half-sarcomere length and the 21 sum of the length from the Z-line to the actin binding site 22 on titin (ℓ 1 ) and the length of the IgD segment that is bound 23 to myosin (L IgD ).Finally, the length of the extra-cellular-24 matrix ℓ ECM is simply half the length of the CE since we are modeling a half-26 sarcomere.

27
We have some freedom to choose the state vector of the 28 model and the differential equations that define how the 29 muscle responds to length and activation changes.The 30 experiments we hope to replicate depend on phenomena 1 that take place at different time-scales: Kirsch et al.'s [5] 2 stochastic perturbations evolve over brief time-scales, while and the elastic-tendon model has five states.For the purpose of comparison, a Hill-type muscle model with a rigid-tendon has a single state (a), while an elastic-tendon model has two states (a and ℓ M ) [17].
Before proceeding, a small note on notation: through- over longer time-scales, where f L (•) is the active-force-16 length curve (Fig. 2A), f PE (•) is the passive-force-length 17 curve (Fig. 2A), and f V (•) is the force-velocity (Fig. 2B).18 The normalized tension developed by the VEXAT model 19 differs from that of a Hill model, Eqn.14, because it has 20 no explicit dependency on ṽ M , includes four passive terms, 21 and a lumped viscoelastic crossbridge element.The four 22 passive terms come from the ECM element f ECM ( l ECM ) 1 (Fig. 3A), the PEVK-IgD element f 2 ( l 2 ) (Fig. 3A and B), the compressive term f KE ( l M ) (prevents l M cos α from reaching a length of 0), and a numerical damping term 4 β ϵ ṽ M (where β ϵ is small).The active force developed by the XE's lumped crossbridge k X o lX + β X o ṽ X is scaled by the 6 fraction of the XE that is active and attached, af L ( lS + LM ), where f L (•) is the active-force-length relation (Fig. 2A).

8
We evaluate f L using lS + LM in Eqn. 15, rather than lM 9 as in Eqn.14, since actin-myosin overlap is independent 10 of crossbridge strain.With f M derived, we can proceed to 11 model the acceleration of the CE, v S , so that it is driven 12 over time by the force imbalance between the XE's active 13 tension and that of a Hill model.

14
We set the first term of v S so that, over time, the CE 15 is driven to develop the same active tension as a Hill-type 16 model [17] (terms highlighted in blue) where τ S is a time constant and f V (v S ) is the force-velocity 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 (f 1 ( l 1 ) and f 2 ( l 2 )) and the ECM (f ECM ( l ECM )).We model titin's elasticity as being due to two serially connected elastic segments: the first elastic segment f 1 ( l 1 ) is formed by lumping together the IgP segment and a fraction Q of the PEVK segment, while the second elastic segment f 2 ( l 2 ) 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 [7] indicate that a Q value of 0.5, positioning the PEVK-actin attachment point that is near the middle of the PEVK segment, allows the model to develop sufficient tension when actively lengthened.The large section of the IgD segment that is bound to myosin is treated as rigid.
The curves that form f ECM ( l ECM ), f 1 ( l 1 ), and f 2 ( l 2 ) have been carefully constructed to satisfy three experimental observations: that the total passive force-length curve of titin and the ECM match the observed passive force-length curve of the muscle (Fig. 2A and Fig. 3A) [58]; that the proportion of the passive force developed by titin and the ECM is within experimental observations [58] (Fig. 3A); and that the Ig domains and PEVK residues show the same relative elongation as observed by Trombitás et al. [26] (Fig. 3C).Even though Trombitás et al.'s [26] measurements come from human soleus titin, we can construct the force-length curves of other titin isoforms if the number of proximal Ig domains, PEVK residues, and distal Ig domains is known (see Appendix ??).Even though Trombitás et al. [26] is from human soleus titin, Since both the passiveforce-length relation and When active muscle is lengthened, it produces an enhanced force that persists long after the lengthening has ceased called residual force enhancement (RFE) [7].For the purposes of the VEXAT model, we assume that RFE is produced by titin.Experiments have demonstrated RFE on both the ascending limb [66] and descending limb of the force-length [7] relation.The amount of RFE depends both 1 on the final length of the stretch [67] and the magnitude of  The passive force-length curve has been decomposed such that 56% of it comes from the ECM while 44% comes from titin to match the average of ECM-titin passive force distribution (which ranges from 43%-76%) reported by Prado et al. [58] (A.).The elasticity of the titin segment has been further decomposed into two serially connected sections: the proximal section consisting of the T12, proximal IgP segment and part of the PEVK segment, and the distal section consisting of the remaining PEVK section and the distal Ig segment (B.).The stiffness of the IgP and PEVK segments has been chosen so that the model can reproduce the movements of IgP/PEVK and PEVK/IgD boundaries that Trombitás et al. [26] (C.) observed in their experiments.That the curves that appear in subplots A. and B. come from scaling the two-segment human soleus titin model to cat soleus.The curves that appear in subplot C compare the human soleus titin model's IgP, PEVK, and IgD force-length relations to the data of Trombitás et al.'s [26] (see in Appendix B for details).the stretch: on the ascending limb the amount of RFE varies 3 with the final length but not with stretch magnitude, while 4 on the descending limb RFE varies with stretch magnitude.

5
To develop RFE, we assume that some point of the PEVK 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, l1 + LT12 , exceeds the length of the actin filament, LA .To break the titin-actin bond at long CE lengths we introduce a smooth step-down function The step-down function u L transitions from one to zero when the titin-actin bond ( l1 + LT12 ) exceeds the actin tip at a length of LA .
The strength of the titin-actin bond also appears to vary nonlinearly with activation.Fukutani and Herzog [68] observed that the absolute RFE magnitude produced by actively lengthened fibers is similar between normal and reduced contractile force states.Since these experiments [68] were performed beyond the optimal CE length, titin 1 could be contributing to the observed RFE as previously 2 described.The consistent pattern of absolute RFE values observed by Fukutani and Herzog [68] could be produced if the titin-actin bond saturated at its maximum strength even at a reduced contractile force state.To saturate the titin-actin bond, we use a final smooth step function where A • is the threshold activation level at which the bond saturates.While we model the strength of the titin-actin bond as being a function of activation, which is proportional Ca 2+ concentration [69], this is a mathematical convenience.The work of Leonard et al. [8] makes it clear that both Ca 2+ and crossbridge cycling are needed to allow titin to develop enhanced forces during active lengthening: no enhanced forces are observed in the presence of Ca 2+ when crossbridge cycling is chemically inhibited.Putting this all together, the active damping acting between the titin and actin filaments is given by the product of u A u S u L β PEVK A , where β PEVK A is the maximum damping coefficient.
With a model of the titin-actin bond derived, we can focus on how the bond location moves in response to applied forces.Since we are ignoring the mass of the titin filament, the PEVK-attachment point is balanced by the forces applied to it and the viscous forces developed between titin and actin due to the active (u A u S u L β PEVK A ) and a small amount of passive damping (β PEVK P ).Since Eqn.20 is linear in ṽ 1 , we can solve directly for it The assumption of whether the tendon is rigid or elastic affects how the state derivative is evaluated and how expensive it is to compute.While all of the position dependent quantities can be evaluated using 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 1 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 elastictendon muscle model is more involved.As is typical of lumped parameter muscle models [16], [17], [70], here we assume that difference in tension, f ϵ , between the CE and 1 the tendon is negligible 10 .During our preliminary simulations it be-3 came clear that treating the tendon as an idealized spring 4 degraded the ability of the model to replicate the experiment 5 of Kirsch et al. [5] particularly at high frequencies.Kirsch 6 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 9 can be accurately reproduced by a spring in parallel with a 10 damper.Due to the way that impedance combines in series 11 11 , the models of both the CE and the tendon need to have 12 parallel spring and damper elements so that the entire MTU, 13 when linearized, appears to be a spring in parallel with a 14 damping element.We model tendon force using a nonlinear 15 spring and damper model where the damping coefficient U k T ( lT ), is a linear scal-17 ing of the normalized tendon stiffness k T by U , a constant 18 scaling coefficient.We have chosen this specific damping 19 model because it fits the data of Netti et al. [72] and cap-20 tures the structural coupling between tendon stiffness and 21 damping (see Appendix B.1 and Fig. 15 for further details).22 Now that all of the terms in Eqn.22 have been explicitly 23 defined, we can use Eqn.22 to solve for v M .Equation 2224 becomes linear in v M after substituting the force models 25 described in Eqns.23 and 15, and the kinematic model 26 equations described in Eqns.8, 9 and 11 (along with the 27 time derivatives of Eqns.[8][9][10][11].After some simplification 28 we arrive at allowing us to evaluate the final state derivative in ẋ.Dur- ing simulation the denominator of ṽ M will always be finite 2 since β ϵ > 0, and α < 90 • due to the compressive element. 3 The evaluation of ẋ in the VEXAT model is free of numeri- cal singularities, giving it an advantage over a conventional 5 Hill-type muscle model [17].In addition, the VEXAT's ẋ  The VEXAT model requires the architectural muscle parameters needed by a conventional Hill-type muscle model as well as additional parameters.The additional parameters are needed for these component models: the lumped viscoelastic XE (Eqn. 1 and 2) , XE-actin dynamics (Eqn.16), the two-segment active titin model (Fig. 3, titin-actin dynamics (Eqns.21), and the tendon damping model (Eqn.23).Fortunately, there is enough experimental data in the literature that values can be found, fitted, or estimated directly for our simulations of experiments on cat soleus (Table 1), and rabbit psoas fibril (see Appendix B for fitting details and Appendix H for rabbit psoas fibril model parameters).

Stochastic Length Perturbation Experiments
In Kirsch et al.'s [5] in-situ experiment, the force response of a cat's soleus muscle under constant stimulation was measured as its length was changed by small amounts.Kirsch et al. [5] applied stochastic length perturbations (Fig. 4A) to elicit force responses from the muscle (in this case a spring-damper Fig. 4B) across a broad range of frequencies (4-90 Hz) and across a range of small length perturbations (1-3.8%ℓ M o ).The resulting time-domain signals can be quite complicated (Fig. 4A and B) but contain rich measurements of how muscle transforms changes in length into changes in force.
As long as muscle can be considered to be linear (a sinsoidal change in length produces a sinsoidal change in force) then system identification methods [78], [79] can be applied to extract a relationship between length x(t) and force y(t) that is hidden in the time domain but clear in the frequency domain.We will give a brief overview of system identification methods here to make methods and results clearer (see Appendix D for details).First, the time-domain signals (x(t) and y(t)) are transformed into an equivalent representation in the frequency-domain (X(s) and Y (s)) as a sum of scaled and shifted sine curves (Fig. 4B and C) using a Fourier transform [78].In the frequency domain, we identify an LTI system of best fit H(s) that describes how muscle transforms changes in length into changes in force such that Y (s) = H(s) X(s).Next, we evaluate how H(s) scales the magnitude (gain) and shifts the timing (phase) of a sinusoid in X(s) into a sinusoid of the same frequency in Y (s) (Fig. 4D).This process is repeated across all frequency-matched pairs of input and output sinusoids to build up a function of how muscle scales (Fig. 4E) and shifts (Fig. 4F) input length sinusoids into output force sinusoids.The resulting transformation turns Table 1: The VEXAT and Hill model's elastic-tendon cat soleus MTU parameters.The VEXAT model uses all of the Hill model's parameters which are highlighted in grey.Short forms are used to indicate: length 'len', velocity 'vel', acceleration 'acc', half 'h', activation 'act', segment 'seg', threshold 'thr', and stiffness 'stiff'.The letters 'R' or 'H' in front of a reference mean the data is from a rabbit or a human, otherwise the data is from cat soleus.The letters following a reference indicate how the data was used to create the parameter: 'C' calculated, 'F' fit, 'E' estimated, or 'S' scaled.Most of the VEXAT model's XE and titin parameters can be used as rough parameter guesses for other muscles because we have expressed these parameters in a normalized space: the values will scale appropriately with changes to ℓ M o and f M o .Titin's force-length curves, however, should be updated if N IgP , N PEVK , or N IgD differ from the values shown below (see Appendix B.3 for details).Finally, the rigid-tendon cat soleus parameters differ from the table below because tendon elasticity affects the fitting of k Each individual input sinsoid is compared with the output sinusoid of the same frequency to evaluate how the system scales and shifts the input to the output (D).This process is repeated across all sinusoid pairs to produce a function that describes how an input sinsoid is scaled (E) and shifted (F) to an output sinusoid using only the measured data (A).
noise and small nonlinearities, the mathematical procedure 7 used to extract H(s) and the corresponding gain and phase profiles is a little more elaborate than we have described 2 (see Appendix D for details).
3 Kirsch et al. [5] used system identification methods to identify LTI mechanical systems that best describes how muscle transforms input length waveforms to output force waveforms.The resulting LTI system, however, is only accurate when the relationship between input and output is approximately linear.Kirsch et al. [5] used the coherence squared, (C xy ) 2 , between the input and output to evaluate the degree of linearity: Y (s) is a linear transformation of X(s) at frequencies in which (C xy ) 2 is near one, while Y (s) cannot be described as a linear function of X(s) at frequencies in which (C xy ) 2 approaches zero.By calculating (C xy ) 2 between the length perturbation and force waveforms, Kirsch et al. [5] identified the bandwidth in which the muscle's response is approximately linear.Kirsch et al. [5] set the lower frequency of this band to 4 Hz, and Fig. 3 of Kirsch et al. [5] suggests that this corresponds to (C xy ) 2 ≥ 0.67 though the threshold for (C xy ) 2 is not reported.The upper frequency of this band was set to the cutoff frequency of the low-pass filter applied to the input (15,35,or 90 Hz).Within this bandwidth, Kirsch et al. [5] compared the response of the specimen to several candidate models and found that a parallel spring-damper best fit the muscle's frequency response.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 varianceaccounted-for (VAF) between the two time-domain signals (25) Astonishingly, Kirsch et al. [5] found that a spring-damper of best fit has a VAF of between 78-99% 12 when compared to the experimentally measured forces f EXP .By repeating this experiment over a variety of stimulation levels (using both electrical stimulation and the crossed-extension reflex) Kirsch et al. [5] showed that these stiffness and damping coefficients vary linearly with the active force developed by the muscle.Further, Kirsch et al. [5] repeated the experiment using perturbations that had a variety of lengths (0.4 mm, 0.8mm, and 1.6mm) and bandwidths (15Hz, 35Hz, and 90Hz) and observed a peculiar quality of muscle: the damping coefficient of best fit increases as the bandwidth of the perturbation decreases (See Figures 3 and 10 of Kirsch et al. [5] for details).Here we simulate Kirsch et al.'s experiment [5] to determine, first, the VAF of the VEXAT model and the Hill model in comparison to a spring-damper of best fit; second, to compare the gain and phase response Figure 5: The perturbation waveforms are constructed by generating a series of pseudo-random numbers, padding the ends with zeros, by filtering the signal using a 2nd order low-pass filter (wave forms with -3dB cut-off frequencies of 90 Hz, 35 Hz and 15 Hz appear in A.) and finally by scaling the range to the desired limit (1.6mm in A.).Although the power spectrum of the resulting signals is highly variable, the filter ensures that the frequencies beyond the -3dB point have less than half their original power (B.). of the models to biological muscle; and finally, to see if 46 the spring-damper coefficients of best fit for both models 1 increase with active force in a manner that is similar to the 2 cat soleus that Kirsch et al. studied [5].

3
To simulate the experiments of Kirsch et al. [5] we begin 4 by creating the 9 stochastic perturbation waveforms used in 5 the experiment that vary in perturbation amplitude (0.4mm 6 ,0.8mm, and 1.6mm) and bandwidth (0-15 Hz, 0-35 Hz, and 7 0-90 Hz) 13 .The waveform is created using a vector that is 8 composed of random numbers with a range of [−1, 1] that 9 begins and ends with a series of zero-valued padding points.Next, a forward pass of a 2 nd order Butterworth filter is 11 applied to the waveform and finally the signal is scaled to 12 the appropriate amplitude (Fig. 5).The muscle model is  varies with active force, we repeated these simulations at 18 ten evenly spaced tensions from 2.5N to 11.5N.Ninety 19 13 For brevity we will refer to the -3 dB frequency of the perturbation waveform rather than the entire bandwidth simulations are required to evaluate the nine different per-20 turbation waveforms at each of the ten tension levels.The time-domain length perturbations and force responses of the modeled muscles are used to evaluate the coherence squared of the signal, gain response, and phase responses of the models in the frequency-domain.Since the response of the models might be more nonlinear than biological muscle, we select a bandwidth that meets (C xy ) 2 > 0.67 but otherwise follows the bandwidths analyzed by Kirsch et al. [5] (see Appendix D for details).
When coupled with an elastic-tendon, the 15 Hz perturbations show that neither model can match the VAF of Kirsch et al.'s analysis [5] (compare Fig. 6A to G), while at 90Hz the VEXAT model reaches a VAF of 92% (Fig. 6D) which is within the range of 78-99% reported by Kirsch et al. [5].In contrast, the Hill model's VAF at 90 Hz remains low at 61% (Fig. 6J).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. 6B to H and E to K), both models have a greater phase shift than Kirch et al.'s data [5] at low frequencies (compare Fig. 6C to I and F to L).The phase response of the VEXAT model to the 90 Hz perturbation (Fig. 6F) shows the consequences of Eqn.16: at low frequencies the phase response of the VEXAT model is similar to that of the Hill model, while at higher frequencies the model's response becomes similar to a spring-damper.This frequency dependent response is a consequence of the first term in Eqn.16: the value of τ S causes the response of the model to be similar to a Hill model at lower frequencies and mimic a spring-damper at higher frequencies.Both models show the same perturbation-dependent phase-response, as the damping coefficient of best fit increases as the perturbation bandwidth decreases: compare the damping coefficient of best fit for the 15Hz and 90Hz profiles for the VEXAT model (listed on   6: The 15 Hz perturbations show that the VEXAT model's performance is mixed: in the time-domain (A.) the VAF is lower than the 78-99% analyzed by Kirsch et al. [5]; the gain response (B.) follows the profile in Figure 3 of Kirsch et al. [5], while the phase response is differs (C.).The response of the VEXAT model to the 90 Hz perturbations is much better: a VAF of 91% is reached in the time-domain (D.), the gain response follows the response of the cat soleus analyzed by Kirsch et al. [5], while the phase-response follows biological muscle closely for frequencies higher than 30 Hz.Although the Hill's time-domain response to the 15 Hz signal has a higher VAF than the VEXAT model (G.), the RMSE of the Hill model's gain response (H.) and phase response (I.) shows it to be more in error than the VEXAT model.While the VEXAT model's response improved in response to the 90 Hz perturbation, the Hill model's response does not: the VAF of the time-domain response remains low (J.), neither the gain (K.) nor phase responses (L.) follow the data of Kirsch et al. [5].Note that the Hill model's 90 Hz response was so nonlinear that the lowest frequency analyzed had to be raised from 4 Hz to 7 Hz to satisfy the criteria that (C xy ) 2 ≥ 0.67. the elastic-tendon Hill model, this improvement has been 51 made by using amounts of damping that exceeds that of biological muscle [5].
The gain and phase profiles of both models deviate from the spring-damper of best fit due to the presence of nonlinearities, even for small perturbations.Some of the VEXAT Figure 7: When coupled with an elastic-tendon, the stiffness (A.) and damping (B.) coefficients of best fit of both the VEXAT model and a Hill model increase with the tension developed by the MTU.However, both the stiffness and damping of the elastic-tendon Hill model are larger than Kirsch et al.'s coefficients (from Figure 12 of [5]), particularly at higher tensions.When coupled with rigid-tendon the stiffness (C.) and damping (D.) coefficients of the VEXAT model remain similar, as the values for k X o and β X o have been chosen to take the tendon model into account.In contrast, the stiffness and damping coefficients of the rigid-tendon Hill model differ dramatically from the elastic-tendon Hill model: while the elastic-tendon Hill model is too stiff and damped, the rigid-tendon Hill model is not stiff enough (compare A. to C.) and far too damped (compare B. to D.).Coupling the Hill model with a rigid-tendon rather than an elastic-tendon does have an unexpected benefit: the VAF increased from the 25%-45% range shown in A. and B. to 87% shown in C. and D. Although the Hill model's time-domain response improves when coupled with a rigid-tendon, the improvement is made with stiffness and damping that deviates from biological muscle [5].also suffers from high degrees of nonlinearity for small 5 perturbations about v M = 0 since the slope of df V dv M is 1 positive and large when shortening, and positive and small 2 when lengthening (Fig. 2B).In contrast, the stiffness and damping of the VEXAT's CE does not change so drastically be- Kirsch et al's data [5] indicates (Fig. 7C and D).VAF in response to the 1.6 mm 35 Hz perturbations (Fig. 8B).It is unclear if biological muscle displays systematic shifts in VAF since Kirsch et al. [5] did not report the VAF of each trial.

Active lengthening on the descending limb
We now turn our attention to the active lengthening in-situ experiments of Herzog and Leonard [7].During these experiments cat soleus was actively lengthened by modest amounts (7-21% ℓ M o ) starting on the descending limb of the active-force-length curve (ℓ M /ℓ M o > 1 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 ℓ M o .Herzog and Leonard's [7] 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 G for details).
When Herzog and Leonard's [7] active ramp-lengthening (Fig. 9A) experiment is simulated, both models produce a force transient initially (Fig. 9B), but for different reasons.The VEXAT model's transient is created when the lumped crossbridge spring (the k X o lX 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.
After the initial force transient, the response of the two models diverges (Fig. 9B): the VEXAT model continues to develop increased tension as it is lengthened, while the Hill model's tension drops before recovering.The VEXAT model's continued increase in force is due to the titin model: when activated, a section of titin's PEVK region remains approximately fixed to the actin element (Fig. 1C).As a result, the ℓ 2 element (composed of part of PEVK segment and the distal Ig segment) continues to stretch and gener-Figure 8: Kirsch et al. [5] noted that the VAF of the spring-damper model of best fit captured between 78-99% across all experiments.We have repeated the perturbation experiments to evaluate the VAF across a range of conditions: two different tendon models, three perturbation bandwidths (15 Hz, 35 Hz, and 90 Hz), three perturbation magnitudes (0.4 mm, 0.8 mm and 1.6 mm), and ten nominal force values (spaced evenly between 2.5N and 11.5N).Each bar in the plot shows the mean VAF across all ten nominal force values, with the whiskers extending to the minimum and maximum value that occurred in each set.The mean VAF of the VEXAT model changes by up to 30% depending on the condition, with the lowest VAF occurring in response to the 15 Hz 1.6 mm perturbation with an elastic tendon (A.), and the highest VAF occurring in response to the 90 Hz perturbations with the rigid tendon (C.).The mean VAF of the Hill model varies by up to 60% depending on the condition, with the lowest VAF values occurring in the 35 Hz 1.6mm trial with the elastic tendon (B.), and the highest VAF values taking place in the 15 Hz 0.4mm trial with the rigid tendon (A.).
ates higher forces than it would if the muscle were being 50 passively stretched.While both the elastic and rigid-tendon versions of the VEXAT model produce the same stereotypical ramp-lengthening response (Fig. 9B), the rigid-tendon model develops slightly more tension because the strain of the MTU is solely borne by the CE.
In contrast, the Hill model develops less force during lengthening when it enters a small region of negative stiffness (Fig. 9B and C) because the passive-force-length curve is too compliant to compensate for the negative slope of the active force-length curve.Similarly, the damping coefficient of the Hill model drops substantially during lengthening (Fig. 9D).Equation 27and 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. 9B), the tension developed by the cat soleus recovers slowly, following a profile that looks strikingly like a first-order decay.The large damping coefficient acting between the titin-actin bond slows the force recovery of the VEXAT model.We have tuned the value of  Once activation is allowed to return to zero, Herzog and 4 Leonard's data shows that the cat soleus continues to de-5 velop a tension that is ∆f B above passive levels (Fig. 9B 6 for t > 8.5s).The force ∆f B is known as passive force 7 enhancement, and is suspected to be caused by titin bind-8 ing to actin [76].Since we model titin-actin forces using One of the great challenges that remains is to decompose 17 how much tension is developed by titin (Fig. 1C) separately 18 from myosin (Fig. 1B) in an active sarcomere.is reached, at which point the titin-actin bond is pulled off 12 the end of the actin filament and the active force is reduced 13 to its passive value.

14
The WLC titin model is not able to reach the extreme 15 lengths observed by Leonard et al. [8].The distal segment 16 of the WLC titin model approaches its contour length early 17 in the simulation and ensures that the the titin-actin bond 18 is dragged off the end of the actin filament at 1.99 ℓ M o 19 (Fig. 10B).After 1.99ℓ M o (Fig. 10B), the tension of the 20 WLC titin model drops to its passive value but continues 21 to increase until the contour lengths of all of the segments 22 of titin are reached at 2.32 ℓ M o .Comparing the response 23 of the linear model to the WLC titin model two things are 24 clear: the linear titin model more faithfully follows the data 25 of Leonard et al. [8], but does so with titin segment lengths 26 that exceed the maximum contour length expected for the 27 isoform of titin in a rabbit myofibril.

28
This simulation has also uncovered a surprising fact: the 29 myofibrils in Leonard et al.'s [8] experiments do not fail 30 at 2.32 ℓ M o , as would be expected by the WLC model of 31 titin, but instead reach much greater lengths (Fig. 2B).32 Physically, it may be possible for a rabbit myofibril to reach 33 these lengths (without exceeding the contour lengths of the 34 proximal Ig, PEVK, and distal Ig segments) if the bond 35 between the distal segment of titin and myosin breaks down.36 This would allow the large Ig segment, that is normally 37 bound to myosin, to uncoil and continue to develop the 38 forces observed by Leonard et al. [8] [9] force-velocity experiment is simulated (A.), the VEXAT model produces a force-velocity profile (blue dots) that approaches zero more rapidly during shortening than the embedded profile f V (•) (red lines).By scaling v M max by 0.95 the VEXAT model (magenta squares) is able to closely follow the force-velocity curve of the Hill model.While the force-velocity curves between the two models are similar, the time-domain force response of the two models differs substantially (B.).The rigid-tendon Hill model exhibits a sharp nonlinear change in force at the beginning (0.1s) and ending (0.21s) of the ramp stretch.compare the results to the force-length and force-velocity 50 curves that are used in the Hill model and in Eqn.16  The VEXAT model produces forces that differ slightly from the f V that is embedded in Eqn.16 while the Hill model reproduces the curve (Fig. 11).The maximum 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 during Kirsch et al.'s experiments [5] (particularly Fig. 6C., and F.).These small differences can be effectively removed by scaling v M max by s V (Fig. 11A has s V = 0.95) to accommodate for the small decrease in force caused by the viscoelastic XE element.[10] force-length experiments were simulated by first passively lengthening the CE, and next by measuring the active force developed by the CE at a series of fixed lengths.Prior to activation, the passive CE was simulated for a brief period of time in a passive state to reduce any history effects due to the active titin element.To be consistent with Gordon et al.'s [10] experiment, we subtracted off the passive force from the active force before producing the active-force-length profile.

Gordon et al.'s
The simulation of Gordon et al.'s [10] experiment shows that the VEXAT model (Fig. 12A, blue dots) produces a force-length profile that is shifted to the right of the Hill model (Fig. 12A, red line) due to the series elasticity introduced by the XE.We can solve for the size of this rightwards shift by noting that Eqn.16 will drive the lS 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

Discussion
A muscle model is defined by the experiments it can replicate and the mechanisms it embodies.We have developed the VEXAT muscle model to replicate the force response of muscle to a wide variety of perturbations [5], [7], [8] while also retaining the ability to reproduce Hill's force-velocity [9] experiment and Gordon et al.'s [10] force-length experiments.The model we have developed uses two mechanisms to capture the force response of muscle over a large variety of time and length scales: first, a viscoelastic crossbridge element that over brief time-scales appears as a springdamper, and at longer time-scales mimics a Hill-model; second, a titin element that is capable of developing active force during large stretches.
The viscoelastic crossbridge and titin elements we have developed introduce a number of assumptions into the model.While there is evidence that the activationdependent stiffness of muscle originates primarily from the stiffness of the attached crossbridges [62], the origins of the activation-dependent damping observed by Kirsch et al. [5] have not yet been established.We assumed that, since the damping observed by Kirsch et al. [5] varies linearly with activation, the damping originates from the attached crossbridges.Whether this damping is intrinsic or is due to some other factor remains to be established.Next, we have also assumed that the force developed by the XE converges to a Hill model [17] given enough time (Eqn.16).A recent experiment of Tomalka et al. [80] suggests the force developed by the XE might decrease during lengthening rather than increasing as is typical of a Hill model [17].If Tomalka et al.'s [80] observations can be replicated, the VEXAT model will need to be adjusted so that the the XE element develops less force during active lengthening while the active-titin element develops more force.Finally, we have assumed that actin-myosin sliding acceleration (due to crossbridge cycling) occurs when there is a force imbalance between the external force applied to the XE and the internal force developed by the XE as shown in Eqn.16.This assumption is a departure from previous models: Hilltype models [16], [17] assume that the tension applied to the muscle instantaneously affects the actin-myosin sliding velocity; Huxley models [11] assume that the actin-myosin sliding velocity directly affects the rate of attachment and detachment of crossbridges.
The active titin model that we have developed makes as-49 sumptions similar to Rode et al. [38] and Schappacher-Tilp   (ℓ M o and f M o ).In contrast, myosin-actin interactions and some titin-actin interactions are temperature-sensitive [82], 38 [83] and may not scale linearly with size.In Sec.related to submaximal contractions: the shift in the peak of the force-length relation [84], and the scaling of the maximum shortening [85].We hope to include these phenomena in a later version of the VEXAT model to more accurately simulate submaximal contractions.
The model we have proposed can replicate phenomena that occur at a wide variety of time and length scales: Kirsch et al.'s experiments [5] which occur over small time and length scales; and the active lengthening experiments of Herzog and Leonard [7] and Leonard et al. [8] which occur over physiological and supra-physiological length scales.In contrast, we have shown in Sec.3.1 to 3.3 that a Hill-type model compares poorly to biological muscle when the same set of experiments are simulated.We expect that a Huxley model [11] is also likely to have difficulty reproducing Kirsch et al.'s experiment [5] because the model lacks an active damping element.Since titin was discovered [23] long after Huxley's model was proposed [11], a Huxley model will be unable to replicate any experiment that is strongly influenced by titin such as Leonard et al.'s experiment [8].
Although there have been several more recent muscle model formulations proposed, none have the properties to simultaneously reproduce the experiments of Kirsch et al. [5], Herzog and Leonard [7], Leonard et al. [8], Hill [9], and Gordon et al. [10].Linearized impedance models [14], [15] can reproduce Kirsch et al.'s experiments [5], but these models lack the nonlinear components needed to reproduce Gordon et al.'s force-length curve [10] and Hill's forcevelocity curve [9].The models of Forcinito et al. [18], and Tahir et al. [41] have a structure that places a contractile element in series with an elastic tendon.While this is a commonly used structure, at high frequencies the lack of damping in the tendon will drive the phase shift between length and force to approach zero.The measurements and model of Kirsch et al. [5], in contrast, indicate that the phase shift between length and force approaches ninety degrees with increasing frequencies.Though the Hill-type models of Haeufle et al. [20] and Günther et al. [22] have viscoelastic tendons, these models have no representation of the viscoelasticity of the CE's attached crossbridges.Similar to the Hill-type muscle model evaluated in this work [17], it is likely that models of Haeufle et al. [20] and Günther et al. [22] will not be able to match the frequency response of biological muscle.De Groote et al. [51], [52] introduced a short-range-stiffness element in parallel to a Hill model to capture the stiffness of biological muscle.While De Groote et al.'s [51], [52]  [38] C. Rode, T. Siebert, and R. Blickhan, "Titin-induced 54 force enhancement and force depression: A 'sticky-  Since the actin filament length varies across species we label it L A .Across rabbits, cats and human skeletal muscle myosin geometry is consistent [75]: a half-myosin is 0.8µm in length with a 0.1µm bare patch in the middle.Thus at full overlap the average point of attachment is 0.45µm from the M-line, or L A −0.45µm from the Z-line at ℓ M o .The lumped stiffness of the actin-myosin load path of a half-sarcomere is the stiffness of three springs in series: a spring representing a L A − 0.45µm length of actin, a spring representing the all attached crossbridges, and a spring representing a 0.45µm section of myosin.
A single half-myosin can connect to the surrounding 2 six actin filaments through 97.9 crossbridges.A 0.800 3 µm half-myosin has a pair crossbridges over 0.700 µm of 4 its length every 14.3nm which amounts to 97.9 per half-5 myosin [12].Although 97.9 crossbridges does not make 6 physical sense, here we will evaluate the stiffness of the 7 CE assuming that fractional crossbridges exist and that 8 attached crossbridges can be perfectly distributed among 9 the 6 available actin filaments: the alternative calculation is 10 more complicated and produces stiffness values that differ 11 only in the 3 rd significant digit.Assuming a duty cycle of 12 20% [81] (values between 5-90% have been reported [86]), 13 at full actin-myosin overlap there will be 19.6 crossbridges 14 attached to the 6 surrounding actin filaments.Assuming 15 that these 19.6 crossbridges are evenly distributed between 16 the 6 actin filaments, each single actin will be attached to 17 3.26 attached crossbridges.Comparing the actin-myosin and titin stiffness ranges (Fig. 14) makes it clear that the stiffness of actin-myosin with 1 attached crossbridge (AM:Low in Fig. 14) is comparable to the highest stiffness values we have estimated for titin (TA:High in Fig. 14).When all 20% of the available crossbridges are attached (AM:High in Fig. 14), the average stiffness of the actin-myosin load path is roughly one order of magnitude stiffer than the highest stiffness values of titin (TA:High in Fig. 14), and two to three orders of magnitude higher than the lowest stiffness titin load path (TP:Low in Fig. 14).Similarly, the maximum XE stiffness and titin stiffness in this work are separated by roughly an order of magnitude: the cat soleus model has a XE stiffness of 47.9f M o /ℓ M o and maximum active titin stiffness of 15 See main ActinMyosinAndTitinStiffness.m in the elife2023 branch of accompanying code respotiory for details. ) does its stiffness (TP:High and TA:High) become comparable to the actin-myosin with a single attached crossbridge (AM:Low).At higher activations and modest lengths, the stiffness of the actin-myosin load path (AM: High) exceeds the stiffness of titin (TP: Low and TA:Low) by between two and three orders of magnitude.At higher activations and longer lengths, the stiffness of the actin-myosin load path (AM: High) exceeds the stiffness of titin by roughly an order of magnitude (TP:High and TA:High).

B.1 Fitting the tendon's stiffness and damping
Similar to previous work [17], we model the force-length relation of the tendon using a quintic Bézier spline (Fig. 15A) that begins at ( lT , f T ) = (1.0,0) (where lT is tendon length normalized by ℓ T s , and f T is tension normalized by f M o ), ends at ( lT , f T ) = (1.0 + e T toe , f T toe ) with a normalized stiffness of k T , and uses the constants f T toe = 2/3 Loeb [74], e T o 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 e T o (Fig. 15B and C).
The similarity in shape is likely not a coincidence.
The nonlinear characteristics (Fig. 15) tendon originates from its microscopic structure.Tendon is composed of many fiber bundles with differing resting lengths [72].Initially the tendon's fiber bundles begin crimped, but gradually stretch as the tendon lengthens, until finally all fiber bundles are stretched and the tendon achieves its maximum stiffness (Fig. 15B) and damping (Fig. 15C) [72].Accordingly, in Eqn.23 we have described the normalized damping of the tendon as being equal to the normalized stiffness of the tendon scaled by a constant U .To estimate U we have used the measurements of Netti et al. [72] (Fig. 15 B and C) and solved a least-squares problem to arrive at a value of U = 0.057.The resulting damping model (Fig. 15C) fits the measurements of Netti et al. [72] closely.

B.2 Fitting the CE's Impedance
We can now calculate the normalized impedance of the XE using the viscoelastic-tendon model we have constructed and Kirsch et al.'s [5] measurements of the impedance of the entire MTU.Since an MTUis structured with a CE in series with a tendon, the compliance of the MTU is given by where k M AT is the stiffness of the CE in the direction of 1 the tendon.We can calculate k M directly by fitting a line 2 to the stiffness vs tension plot that appears in Figure 12 3 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. 5 Here we use a nominal tension of 5N so that we can later 6 compare our model to the 5N frequency response reported in Figure 3 of Kirsch et al. [5].Since Kirsch et al. [5] we can use a similar procedure to evaluate β M AT , the damp-19 ing of the CE along the tendon.The value of β M that 20 best fits the damping vs. tension plot that appears in Fig- 21 ure 12 of Kirsch et al. [5] at a nominal tension of 5N is 22 0.0198 Ns/mm.The tendon damping model we have just 23 constructed develops 0.697 Ns/mm at a nominal load of 5N.24 Using Eqn.34, we arrive at β M AT =0.020 Ns/mm.Due to 25 the pennation model, the stiffness and damping values of 26 the CE differ from those along the tendon.

27
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 Figure 15: The normalized tendon force-length curve (A) has been been fit to match the cat soleus tendon stiffness measurements of Scott and Loeb [74].The data of Netti et al. [72] allow us to develop a model of tendon damping as a linear function of tendon stiffness.By normalizing the measurements of Netti et al. [72] by the maximum storage modulus we obtain curves that are equivalent to the normalized stiffness (B) and damping (C) of an Achilles tendon from a rabbit.Both normalized tendon stiffness and damping follow similar curves, but at different scales, allowing us to model tendon damping as a linear function of tendon stiffness (C).
which simplifies to Similarly, the constant thickness pennation model means which leads to Recognizing that The nonlinear force-length curves used to describe titin (f 1 ( l 1 ) and f 2 ( l 2 ) in series), and the ECM (f ECM ( l 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 [7]; the proportion of titin's contribution relative to the ECM must be within measured bounds [58]; and finally the stiffness of the f 2 ( l 2 ) must be a linear scaling of f 1 ( l 1 ) to match the observations of Trombitás et al. [26].
First, we fit the passive force-length curve to the data of Herzog and Leonard [7] to serve as a reference.The curve f PE begins at the normalized length and force coordinates of (1+e PE 0 ,0) with a slope of 0, ends at (1+e the squared differences between f PE and the passive forcelength data of Herzog and Leonard [7] (Fig. 2A shows both the data and the fitted f PE curve) are minimized.While f PE is not used directly in the model, it serves as a useful reference for constructing the ECM and titin force-length curves.We assume that the ECM force-length curve is a where P is a constant.In this work, we set P to 56% which is the average ECM contribution that Prado et al. [58] measured across 5 different rabbit skeletal muscles 16 .The remaining fraction, 1 − P, of the force-length curve comes from titin.
In mammalian skeletal muscle, titin has three elastic segments [58] connected in series: the proximal Ig segment, the PEVK segment, and the distal Ig segment that is between the PEVK region and the myosin filament (Fig. 1A).Trombitás et al. [26]  Ig segment that overlaps with myosin is rigid (Fig. 1A), and 43 that during passive stretching the tension throughout the 1 titin filament is uniform.Since the sarcomeres in Trombitás 2 et al.'s [26] experiments were passively stretched it is rea-3 sonable to assume that tension throughout the free part of 4 the titin filament is uniform because the bond between titin 5 and actin depends on calcium [31], [36] and crossbridge 6 attachment [8].

7
We begin by digitizing the data of Figure 5 of Trombitás 8 et al. [26] and using the least-squares method to fit lines to 9 Z ℓ IgP/PEVK and Z ℓ PEVK/IgD (where the superscripts 10 mean from ℓ to and so Z ℓ IgP/PEVK is the distance from the 11 Z-line to the border of the IgP/PEVK segments).From 12 these lines of best fit we can evaluate the normalized length 13 of the proximal Ig segment the normalized length of the PEVK segment and the normalized length of the distal Ig segment given the lM .The coefficients that best fit the data from 19 Trombitás et al. [26] appear in Table 2.  [26] results to a rabbit psoas, in principle this approach can be applied to any isoform of titin provided that its geometry is known, and the Ig domains and PEVK residues in the target titin behave similarly to those in human soleus titin.
The only detail that remains is to establish the shape of the IgP, PEVK, and IgD force-length curves.Studies of individual titin filaments, and of its segments, make it clear that titin is mechanically complex.For example, the tandem Ig segments (the IgD and IgP segments) are composed of many folded domains (titin from human soleus has two Ig segments that together have nearly 100 domains [26]).
Each domain appears to be a simple nonlinear spring until it unfolds and elongates by nearly 25 nm in the process [88].Unfolding events appear to happen individually during lengthening experiments [88], with each unfolding event occurring at a slightly greater tension than the last, giving an Ig segment a force-length curve that is saw-toothed.
Although detailed models of titin exist that can simulate the folding and unfolding of individual Ig domains, this level under passive lengthening.These coefficients have been extracted from data of Figure 5 of Trombitás et al. [26] using a least-squares fit.Since Figure 5 of Trombitás et al. [26] plots the change in segment length of a single titin filament against the change in length of the entire sarcomere, the resulting slopes are in length normalized units.The slopes sum to 0.5, by construction, to reflect the fact that these three segments of titin stretch at half the rate of the entire sarcomere (assuming symmetry).The cat soleus titin segment coefficients have been formed using a simple scaling of the human soleus titin segment coefficients, and so, are similar.Rabbit psoas titin geometry [58] differs dramatically from human soleus titin [26] and produce a correspondingly large difference in the coefficients that describe the length of the segments of rabbit psoas titin.

Coefficient
Human Feline Rabbit Soleus [26]  average sarcomere length at failure 3.38 ℓ M o was longer than 21 the maximum length of 2.4-2.7 ℓ M o that would be predicted from the geometry of rabbit psoas titin 18 .
Since we are interested in a computationally efficient model that is accurate at the whole muscle level, we model titin as a multi-segmented nonlinear spring but omit the states needed to simulate the folding and unfolding of Ig domains.Simulations of active lengthening using our titin model will exhibit the enhanced force development that appears in experiments [7], [8], but will lack the nonlinear saw-tooth force-length profile that is measured when individual titin filaments are lengthened [88].To have the broadest possible application, we will fit titin's force-length curves to provide reasonable results for both moderate [7] and large active stretches [8].Depending on the application, it may be more appropriate to use a stiffer force-length curve for the Ig segment if normalized sarcomere lengths stays within 1.5 ℓ M o and no unfolding events occur as was done by Trombitás et al. [65].
To ensure that the serially connected force-length curves of f 1 ( l 1 ) and f 2 ( l 2 ) closely reproduce (1 − P)f PE ( lM ), we are going to use affine transformations of f PE to describe f 1 ( l 1 ) and f 2 ( l 2 ).The total stiffness of the halfsarcomere titin model is given by which is formed by the series connection of f 1 ( l 1 ) and Since each of titin's segments is exposed to the same tension in Trombitás et al.'s experiment [26] the slopes of the lines that Eqns.53-55 describe are directly proportional to the relative compliance (inverse of stiffness) between of each of titin's segments.Using this fact, we can define the normalized stretch rates of the proximal titin segment 61) 18 Rabbit psoas titin [58] attaches at the Z-line with a 100nm rigid segment that spans to T12 epitope, is followed by 50 Ig domains, 800 PEVK residues, and another 22 Ig domains until it attaches to the 800 nm half-myosin filament which can also be considered rigid.If the Ig domains were all unfolded (adding around 25 nm [88]) and each PEVK residue could reach a maximum length of between 0.32nm [65] (see Fig. 5: 700nm/2174 residues is 0.32 nm per residue) to 0.38 nm [90] (see pg. 254), two titins in series would reach a length of 2(100nm + 72(25nm) + 800(0.32nm-0.38nm)+ 800 nm) = 5192-6008nm.Since rabbit sarcomeres have an ℓ M o of 2.2µm a sarcomere could be stretched to a length between 5192-6008nm, or 2.4-2.7 ℓ M o , before the contour lengths of the tandem Ig and PEVK segments is reached.and the distal titin segment 62) which are proportional to the compliance of two titin seg-1 ments in our model.When both the f 1 ( l 1 ) and f 2 ( l 2 ) 2 curves are beyond the toe region the stiffness is a constant and Dividing Eqn.63 by 64 eliminates the unknown ∆ f and 6 results in an expression that relates the ratio of the terminal 7 linear stiffness of f 1 ( l 1 ) and f 2 ( l 2 ) to the relative changes in Eqns.61 and 62. Substituting Eqns.65, and 60 into Eqn.59 yields the expression which can be simplified to and this expression can be evaluated using the terminal 12 stiffness of titin k Ti * and the coefficients listed in Table 2. 13 Substituting Eqn.67 into Eqn.65 yields The curves f 1 ( l 1 ) and f 2 ( l 2 ) can now be formed by 15 scaling and shifting the total force-length curve of titin (1 − 16 P)f PE .By construction, titin's force-length curve develops 17 a tension of (1 − P), and has reached its terminal stiffness, 18 when the CE reaches a length lM * such that f PE ( lM * ) = 1.By construction, the spring network formed by the f ECM ( l ECM ), f 1 ( l 1 ), and f 2 ( l 2 ) 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 l1 is effectively fixed in place, the distal segment of titin contributes higher forces since l2 undergoes higher strains (Fig. 3A).Finally, when the experiment of Trombitás et al. [26] are simulated the movements of the IgP/PEVK and PEVK/IgD boundaries in the titin model closely follow the data (Fig. 3C).
The curve is linearly extrapolated, so too are the f ECM ( l ECM ), f 1 ( l 1 ), and f 2 ( l 2 ) 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 LW S (the superscript W means WLC) so that the WLC model can made to be continuous with f 1 ( l 1 ) and f 2 ( l 2 ).The normalized force developed by our WLC model is given by where B is a scaling factor and the normalized segment length lW is defined as where LW S is the slack length, and L W C is the contour length of the segment.To extend the f 1 ( l Next, we define the slack length by linearly extrapolating 8 backwards from the final fitted force (1 − P) and similarly We can now solve for B in Eqn.71 so that f 1 ( l 1 ) and 11 f 2 ( l 2 ) are continuous with each respective WLC extrap-12 olation.However, we do not use the WLC model directly 13 because it contains a numerical singularity which is prob-14 lematic during numerical simulation.Instead, we add an 15 additional Bézier segment to fit the WLC extension that 16 spans between forces of (1 − P) and twice the normal-17 ized failure force (2 × 5.14f M o ) noted by Leonard et al. 18 [8].To fit the shape of the final Bézier segment, we adjust 19 the locations of the internal control points to minimize the 20 squared differences between the modified WLC model and 21 the final Bézier curve (Fig. 10A).The final result is a set 22 of curves (f 1 ( l 1 ), f 2 ( l 2 ), and f ECM ( l ECM )) which, be-23 tween forces 0 and (1 − P), will reproduce f PE , Trombitás 24 et al.'s measurements [26], and do so with a reasonable 25 titin-ECM balance [58].For forces beyond (1 − P), the 26 curve will follow the segment-specific WLC model up to 27 twice the expected failure tension noted by Leonard et By noting that all of our chosen state variables in Eqn.
As with k M , the expression for β M can be reduced to

C Model Initialization
Solving for an initial state is challenging since we are given a, ℓ P , and v P and must solve for v S , ℓ S , and ℓ 1 for a rigidtendon 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 35 well approximated as constant.Our preliminary simulations found that satisfactory solutions were found by iterating over both lM and ṽ M using a nested bisection search that looks for values which approximately satisfies Eqn.22, result in small values for v S 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 lM .Given a, ℓ P , v P , and a candidate value of lM , we can immediately solve for α and ℓ T using the pennation model.We can numerically solve for the value of another state, ℓ 1 , using the kinematic relationship between ℓ M and ℓ 1 and by assuming that the two titin segments are in a force equilibrium In the inner loop, we iterate over values of ṽ M 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 ṽ M that best satisfies Eqn.
22. Prior to evaluating Eqn.22, we need to set both ṽ X and lX .Here we choose a value for ṽ X that will ensure that the XE is not producing transient forces and we use fixed-point iteration to solve for lX such that Eqn.16 evaluates to zero.Now the value of ṽ S can be directly evaluated using the candidate value of ṽ M , the first derivative of Eqn. 9, and the fact that we have set ṽ X to zero.
Finally, the error of this specific combination of lM and ṽ M is evaluated using Eqn.22, where the best solution leads to the lowest absolute value for of f ϵ in Eqn.22.If a rigidtendon model is being initialized the procedure is simpler because the inner loop iterating over ṽ M is unnecessary: given v P and ṽ X are zero, the velocities ṽ M and ṽ S 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 ( lMT ) and nominal force where * is the convolution operator.Each of these signals In Eqn.87  24 Though Koopmans's [79] estimator is a great improve-25 ment over Eqn.88, the accuracy of the estimate can be fur-26 ther improved using Welch's method [91].Welch's method 27 [91] breaks up the time domain signal into K segments, 28 transforms each segment into the frequency domain, and 29 returns the average across all segments.Using Welch's 30 method [91]  F Supplementary plots: Gain and phase response rigid-tendon muscle models   3) [5] than when an elastic-tendon is used.This improvement in accuracy is also observed at the 90 Hz perturbation (D., E., and F.), though the phase response of the model departs from Kirsch et al.'s data [5] for frequencies lower than 30 Hz. Parts of the Hill model's response to the 15 Hz perturbation are better with a rigid-tendon, with a higher VAF (G.), a lower RMSE gain-response (H.).but have a poor phase-response (I.).In response to the higher frequency perturbations, the Hill model's response is poor with an elastic (see Fig. 6) or rigid-tendon.The VAF in response to the 90 Hz perturbation remains low (J.), and neither the gain (K.) nor the phase response of the Hill model (L.) follow the data of Kirsch et al. [5].The rigid-tendon Hill model's nonlinearity was so strong that the lowest frequency analyzed had to be raised from 4 Hz to 21 Hz to meet the criteria that (C xy ) 2 ≥ 0.67.As with the prior simulations the Hill model exhibits a small region of negative stiffness introduced by the descending limb of the force-length curve (C.) and a drop in damping (D.).Note: neither model was fitted to this trial.

Figure 1 :
Figure 1: The name of the VEXAT model comes from the viscoelastic crossbridge and active titin elements (A.) in the model.Active tension generated by the lumped crossbridge flows through actin, myosin, and the adjacent sarcomeres to the attached tendon (B.).Titin is modeled as two springs of length ℓ 1 and ℓ 2 in series with the rigid segments L T12 and L IgD .Viscous forces act between the titin and actin in proportion to the activation of the muscle (C.), which reduces to negligible values in a purely passive muscle (D.).We modeled actin and myosin as rigid elements; the XE, titin, and the tendon as viscoelastic elements; and the ECM as an elastic element.

37
the VEXAT model due to the viscoelastic (VE) crossbridge 38 (X) and active-titin (AT) elements of the model.39 To reduce the number of states needed to simulate the 40 VEXAT model, we lump all of the attached crossbridges 41 into a single lumped crossbridge element (XE) that attaches 42 at ℓ S (Fig. 1A) and has intrinsic stiffness and damping 43 properties that vary with the activation and force-length 44 properties of muscle.The active force developed by the XE 45 at the attachment point to actin is transmitted to the main 1 myosin filament, the M-line, and ultimately to the tendon 2 (Fig. 1B).In addition, since the stiffness of actin [60] and 3 4 myosin filaments [61] greatly exceeds that of crossbridges [62], we treat actin and myosin filaments as rigid to reduce the number of states needed to simulate this model.
linearly with f M o and inversely with ℓ M o .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. [58].

Figure 2 :
Figure2: The model relies on Bézier curves to model the nonlinear effects of the active-force-length curve, the passiveforce-length curves (A.), and the force-velocity curve (B.).Since nearly all of the reference experiments used in Sec. 3 have used cat soleus, we have fit the active-force-length curve (f L (•)) and passive-force-length curves (f PE (•)) to the cat soleus data of Herzog and Leonard 2002[7].The concentric side of the force-velocity curve (f V (•)) has been fitted to the cat soleus data ofHerzog and Leonard 1997 [63].

2 8 For
readers who require an activation model with continuity to the second-derivative, the model of De Groote et al. [64] is recommended. 2

Figure 3 :
Figure3: The passive force-length curve has been decomposed such that 56% of it comes from the ECM while 44% comes from titin to match the average of ECM-titin passive force distribution (which ranges from 43%-76%) reported by Prado et al. [58] (A.).The elasticity of the titin segment has been further decomposed into two serially connected sections: the proximal section consisting of the T12, proximal IgP segment and part of the PEVK segment, and the distal section consisting of the remaining PEVK section and the distal Ig segment (B.).The stiffness of the IgP and PEVK segments has been chosen so that the model can reproduce the movements of IgP/PEVK and PEVK/IgD boundaries that Trombitás et al.[26] (C.) observed in their experiments.That the curves that appear in subplots A. and B. come from scaling the two-segment human soleus titin model to cat soleus.The curves that appear in subplot C compare the human soleus titin model's IgP, PEVK, and IgD force-length relations to the data of Trombitás et al.'s[26] (see in Appendix B for details).

6
does not require iteration to numerically solve a root, giving 7 it an advantage over a singularity-free formulation of the 8 Hill model [17].As with previous models, initializing the 9 model's state is not trivial and required the derivation of a 10 model-specific method (see Appendix C for details).

11 3 Biological Benchmark Simulations 12 In 15 3. 24 task becomes familiar. The active lengthening experiment 25 of Herzog and Leonard [ 7 ] shows that modestly stretched 26 ( 7 - 29 [ 17 ]
order to evaluate the model, we have selected three ex-13 periments that capture the responses of active muscle to 14 small, medium, and large length changes.The small (1-8% ℓ M o ) stochastic perturbation experiment of Kirsch et al.[5] demonstrates that the impedance of muscle is well 17 described by a stiff spring in parallel with a damper, and 18 that the spring-damper coefficients vary linearly with active 19 force.The active impedance of muscle is such a fundamen-20 tal part of motor learning that the amount of impedance, 21 as indicated by co-contraction, is used to define how much 22 learning has actually taken place [1], [77]: co-contraction 23 is high during initial learning, but decreases over time as a 21% ℓ M o ) biological muscle has positive stiffness even 27 on the descending limb of the active force-length curve 28 (ℓ M o > 1).In contrast, a conventional Hill model [16], can have negative stiffness on the descending limb 30 of the active-force-length curve, a property that is both 31 mechanically unstable and unrealistic.The final active 32 lengthening experiment of Leonard et al. [8] unequivo-33 cally demonstrates that the CE continues to develop ac-34 tive forces during extreme lengthening (329% ℓ M o ) which 35 exceeds actin-myosin overlap.Active force development 36 beyond actin-myosin overlap is made possible by titin, and 37 its activation dependent interaction with actin [8].The 38 biological benchmark simulations conclude with a replication of the force-velocity experiments of Hill [9] and the force-length experiments of Gordon et al. [10].

4 35 6 (Figure 4 :
Figure 4: Evaluating a system's gain and phase response begins by applying a pseudo-random input signal to the system and measuring its output (A).Both the input and output signals (A) are transformed into the frequency domain by expressing these signals as an equivalent sum of scaled and shifted sinusoids (simple example shown in B and C).Each individual input sinsoid is compared with the output sinusoid of the same frequency to evaluate how the system scales and shifts the input to the output (D).This process is repeated across all sinusoid pairs to produce a function that describes how an input sinsoid is scaled (E) and shifted (F) to an output sinusoid using only the measured data (A).

13 then
activated until it develops a constant tension at a length 14 of ℓ M o .The musculotendon unit is then simulated as the 15 length is varied using the previously constructed waveforms 16 while activation is held constant.To see how impedance

17
Fig. 6A.and D.) and the Hill model (listed on Fig. 6G.and J., respectively).The closeness of each model's response to the springdamper of best fit changes when a rigid-tendon is used instead of an elastic-tendon.While the VEXAT model's response to the 15 Hz and 90 Hz perturbations improves slightly (compare to Fig. 6A-F to Fig. 16A-F in Appendix F), the response of the Hill model to the 15 Hz perturbation changes dramatically with the time-domain VAF rising from 51% to 86% (compare to Fig. 6G-L to Fig. 16G-L in Appendix F).Although the Hill model's VAF in response to the 15 Hz perturbation improved, the frequency response contains mixed results: the rigid-tendon Hill model's gain response is better (Fig. 16H in Appendix F), while the phase response is worse in comparison to the elastic-tendon Hill model.While the rigid-tendon Hill model produces a better time-domain response to the 15 Hz perturbation than

Figure
Figure6: The 15 Hz perturbations show that the VEXAT model's performance is mixed: in the time-domain (A.) the VAF is lower than the 78-99% analyzed by Kirsch et al.[5]; the gain response (B.) follows the profile in Figure3of Kirsch et al.[5], while the phase response is differs (C.).The response of the VEXAT model to the 90 Hz perturbations is much better: a VAF of 91% is reached in the time-domain (D.), the gain response follows the response of the cat soleus analyzed by Kirsch et al.[5], while the phase-response follows biological muscle closely for frequencies higher than 30 Hz.Although the Hill's time-domain response to the 15 Hz signal has a higher VAF than the VEXAT model (G.), the RMSE of the Hill model's gain response (H.) and phase response (I.) shows it to be more in error than the VEXAT model.While the VEXAT model's response improved in response to the 90 Hz perturbation, the Hill model's response does not: the VAF of the time-domain response remains low (J.), neither the gain (K.) nor phase responses (L.) follow the data of Kirsch et al.[5].Note that the Hill model's 90 Hz response was so nonlinear that the lowest frequency analyzed had to be raised from 4 Hz to 7 Hz to satisfy the criteria that (C xy ) 2 ≥ 0.67.
model's nonlinearities in this experiment come from the tendon model (compare to Fig. 6A-F to Fig. 16A-F in Appendix F), since the response of the VEXAT model with a rigid-tendon stays closer to the spring-damper of best fit.The Hill model's nonlinearities originate from the underlying expressions for stiffness and damping of the Hill model, which are particularly nonlinear with a rigid tendon model (Fig. 16G-L in Appendix F) The stiffness of a Hill model's

1 19 ( 23 (
cause neither of these terms depend on the slope of the 2 force-length relation, nor the force-velocity relation (see across a range of isometric forces, Kirsch et al.[5] were 6 able to show that the stiffness and damping of a muscle 7 varies linearly with the active tension it develops (see Fig-8ure 12 of[5]).We have repeated our simulations of Kirsch 9 et al.'s [5] experiments at ten nominal forces (spaced evenly 10 between 2.5N and 11.5 N) and compared how the VEXAT 11 model and the Hill model's stiffness and damping coeffi-12 cients compare to Figure 12 of Kirsch et al. [5] (Fig. 7).The 13 VEXAT model develops stiffness and damping similar to 14 Kirsch et al.'s data when coupled with either a viscoelastic-15 tendon (Fig. 7A & B) or a rigid-tendon (Fig. 7C & D), and 16 with a high VAF.In contrast, when the Hill model is cou-17 pled with an elastic-tendon both its stiffness and damping 18 are larger than Kirsch et al.'s data [5] at the higher tensions Fig. 7A and B).This pattern changes when simulating a 20 Hill model with a rigid-tendon: the model's stiffness is now 21 zero (Fig. 7C), while the model's final damping coefficient 22 is nearly three times the value measured by Kirsch et al.Fig. 7D).Though a stiffness of zero may seem surprising, 24 Eqn.26 shows a stiffness of zero is possible at ℓ M o the 25 nominal CE length of these simulations: at ℓ M o the slope 26 of the active force-length curve is zero and the slope of the 27 passive force-length curve is negligible.The tendon model 28 also affects the VAF of the Hill model to a large degree: the 29 elastic-tendon Hill model has a low VAF 25%-45% (Fig. 30 7A & B) while the rigid-tendon Hill model has a much 31 higher VAF of 87%.Although the VAF of the rigid-tendon 32 Hill model is acceptable, these forces are being generated 33 in a completely different manner than biological muscle as 34

35 39 (
When the VAF of the VEXAT and Hill model is evalu-36 ated across a range of nominal tensions (ten values from 37 2.5 to 11.5N), frequencies (15 Hz, 35 Hz, and 90 Hz), am-38 plitudes (0.4mm, 0.8mm, and 1.6mm), and tendon types rigid and elastic) two things are clear: first, that the VAF 40 of 79-99% of the VEXAT model is closer to the range of 41 78-99% reported by Kirsch et al. [5] than the range of 25-42 87% of the Hill model (Fig. 8); and second, that there are 43 systematic variations in VAF, stiffness, and damping across 44 the different perturbation magnitudes and frequencies (see Tables 5 and 5 in Appendix E).Both models produce worse 46 VAF values when coupled with an elastic-tendon (Fig. 8A, 47 B, and C), though the Hill model is affected most: the VAF 48of the elastic tendon Hill model is 60% lower than the VAF of the rigid tendon model for low frequencies perturbations (Fig.8A and B).While the VEXAT model's lowest VAF occurs in response to the low frequency perturbations (Fig.8A), the Hill model's lowest VAF varies with both tendon type and frequency: the rigid-tendon Hill model has its lowest VAF in response to the 1.6 mm 90 Hz perturbations (Fig.8C) while the elastic-tendon Hill model has its lowest 8f M o /(ℓ M o /s) for the elastic-tendon model, and 77.7f M o /(ℓ M o /s) for the rigid-tendon model, to match the rate of force decay of the cat soleus in Herzog and Leonard's data [7].The Hill model, in contrast, recovers to its isometric value quite rapidly.Since the Hill model's 27 force enhancement during lengthening is a function of the 1 rate of lengthening, when the lengthening ceases, so too 2 does the force enhancement.

Figure 9 :
Figure9: Herzog and Leonard[7] actively lengthened (A.) cat soleus on the descending limb of the force-length curve (where l M > 1 in Fig.2A) and measured the force response of the MTU (B.).After the initial transient at 2.4s the Hill model's output force drops (B.) because of the small region of negative stiffness (C.) created by the force-length curve.In contrast, the VEXAT model develops steadily increasing forces between 2.4 − 3.4s and has a consistent level of stiffness (C.) and damping (D.).

Figure 10 : 71 . 9 2 4 VEXAT model with the linear titin element can match the 5 experimental 6 model cannot produce active force for lengths greater than 7 1
Figure 10: In the VEXAT model we consider two different versions of the force-length relation for each titin segment (A): a linear extrapolation, and a WLC model extrapolation.Leonard et al. [8] observed that active myofibrils continue to develop increasing amounts of tension beyond actin-myosin overlap (B, grey lines with ±1 standard deviation shown).When this experiment is replicated using the VEXAT model (B., blue & magenta lines) and a Hill model (C.red lines), only the VEXAT model with the linear extrapolated titin model is able to replicate the experiment with the titin-actin bond slipping off of the actin filament at 3.38ℓ M o .

Figure 11 :
Figure 11: When Hill's[9] force-velocity experiment is simulated (A.), the VEXAT model produces a force-velocity profile (blue dots) that approaches zero more rapidly during shortening than the embedded profile f V (•) (red lines).By scaling v M max by 0.95 the VEXAT model (magenta squares) is able to closely follow the force-velocity curve of the Hill model.While the force-velocity curves between the two models are similar, the time-domain force response of the two models differs substantially (B.).The rigid-tendon Hill model exhibits a sharp nonlinear change in force at the beginning (0.1s) and ending (0.21s) of the ramp stretch.

Figure 12 :X o to 3 the
Figure 12: When Gordon et al.'s [10] passive and active force-length experiments are simulated, the VEXAT model (blue dots) and the Hill model (red lines) produce slightly different force-length curves (A.) and force responses in the time-domain (B.).The VEXAT model produces a right shifted active force-length curve, when compared to the Hill model due to the series elasticity of the XE element.By shifting the underlying curve by 2 k X o to the left the VEXAT model (magenta squares) can be made to exactly match the force-length characteristic of the Hill model.

1 8 [ 39 ] 14 [ 5 ]
et al [40]: some part of the PEVK segment bonds actin, and 2 this bond cannot do any positive work on titin.The assump-3 tion that the bond between titin and actin cannot do positive 4 work means that titin cannot be actively preloaded: it can 5 only develop force when it is stretched.In contrast, two 6 mechanisms have been proposed that make it possible for 7 titin to be preloaded by crossbridge cycling: Nishikawa's winding filament theory and DuVall et al.'s [44] titin en-9 tanglement hypothesis.If titin were significantly preloaded 10 by crossbridge cycling, the titin load path would support 11 higher forces and the myosin-actin load path would bear 12 less force.Accordingly, the overall stiffness of the CE 13 would be reduced, affecting our simulations of Kirsch et al. : lower myosin-actin loads mean fewer attached cross-15 bridges, since crossbridges are stiff in comparison to titin 16 the stiffness of the CE would fall (see Appendix A).Hope-17 fully experimental work will clarify if titin can be actively 18 preloaded by crossbridges in the future.19 Both the viscoelastic crossbridge and active titin elements 20 include simple myosin-actin and titin-actin bond models 21 that improve accuracy but have limitations.First, the vis-22 coelastic crossbridge element has been made to represent 23 a population of crossbridges in which the contribution of 24 any single crossbridge is negligible.Though it may be pos-25 sible for the XE model to accurately simulate a maximally 26 activated single sarcomere (which has roughly 20 attached 27 crossbridges per half sarcomere [12], [81]) the accuracy 28 of the model will degrade as the number of attached cross-29 bridges decreases.When only a single crossbridge remains 30 the XE model's output will be inaccurate because it can 31 only generate force continuously while a real crossbridge 32 generates force discretely each time it attaches to, and de-33 taches from, actin.Next, we have used two equations, Eqns.

34 16 and 21 ,
that assume myosin-actin and titin-actin interac-35 tions are temperature-invariant and scale linearly with size 36

3 [ 26 (
formulation improves upon a Hill model it is unlikely to reproduce Kirsch et al.'s experiment[5] because we have shown in Sec.3.1 that a Hill model has a frequency response that differs from biological muscle.Rode et al.'s [38] muscle model also 51 uses a Hill model for the CE and so we expect that this 1 model will have the same difficulties reproducing Kirsch 2 et al.'s [5] experiment.Schappacher-Tilp et al.'s model 40] extends a Huxley model [11] by adding a detailed4 titin element.Similar to a Huxley model, Schappacher-Tilp 5 et al.'s model [40] will likely have difficulty reproducing 6 Kirsch et al.'s experiment [5] because it is missing an active 7 damping element.8 While developing this model, we have come across open 9 questions that we hope can be addressed in the future.How 10 do muscle stiffness and damping change across the force-11 length curve?Does stiffness and damping change with 12 velocity?What are the physical origins of the active damp-13 ing observed by Kirsch et al. [5]?What are the conditions 14 that affect passive-force enhancement, and its release?In 15 addition to pursuing these questions, we hope that other 16 researchers continue to contribute experiments that are 17 amenable to simulation, and to develop musculotendon 18 models that overcome the limitations of our model.To help 19 others build upon our work, we have made the source code 20 of the model and all simulations presented in this paper 21 available online 14 .Financial support is gratefully acknowledged from the 24 Deutsche Forschungsgemeinschaft (DFG, German Re-25 search Foundation) under Germany's Excellence Strategy EXC 2075 -390740016) through the Stuttgart Center for 27 Simulation Science (SimTech), from DFG grant no.MI 28 2109/1-1, the Lighthouse Initiative Geriatronics by StMWi 29 Bayern (Project X, grant no.5140951), and NSERC of 30 Canada.

Figure 13 :
Figure13: To evaluate the stiffness of the actin-myosin load path, we first determine the average point of attachment.Since the actin filament length varies across species we label it L A .Across rabbits, cats and human skeletal muscle myosin geometry is consistent[75]: a half-myosin is 0.8µm in length with a 0.1µm bare patch in the middle.Thus at full overlap the average point of attachment is 0.45µm from the M-line, or L A −0.45µm from the Z-line at ℓ M o .The lumped stiffness of the actin-myosin load path of a half-sarcomere is the stiffness of three springs in series: a spring representing a L A − 0.45µm length of actin, a spring representing the all attached crossbridges, and a spring representing a 0.45µm section of myosin.

18
At full overlap, the Z-line is 1 actin filament length L A 19 (1.12 µm in rabbits[60]) from the M-line.The average 20 point of crossbridge attachment is in the middle of the half-21 myosin at a distance of 0.45 µm from the M-line (0.1 µm is 22bare and 0.35 µm is half of the remaining length), which is 23 L A −0.45 µm from the Z-line.A single actin filament has a stiffness of 46 − 68 pN/nm [60] while a single crossbridge has a stiffness of 0.69 ± 0.47 pN/nm [62].Since stiffness scales inversely with length, actin's stiffness between the Zline and the average attachment point is 81.8−121 pN/nm.Finally, the stiffness of each actin filament and its 3.26 attached crossbridges is 0.712 − 3.67 pN/nm and all 6 together have a stiffness of 4.27 − 22.0 pN/nm.Myosin has a similar stiffness as a single actin filament [61], with the section between the average attachment point and the M-line having a stiffness of 76.9 − 113 pN/nm.The final active stiffness of half-sarcomere is 4.05 − 18.4 pN/nm which comes from comes from the series connection of the group of 6 actin filaments, with 19.6 crossbridges, and finally the single myosin filament.When this procedure is repeated assuming that only a single crossbridge is attached the stiffness drops to 0.22−1.15pN/nm, which is slightly less than the stiffness of a single crossbridge 15 .The force-length profile of a single rabbit titin has been measured by Kellermayer et al. [87] using laser tweezers to apply cyclical stretches.By digitizing Fig. 4B (blue line) of Kellermayer et al. [87] we arrive at a stiffness for titin of 0.0058 − 0.0288 pN/nm at 2 µm (for a total sarcomere length of 4 µm or 1.62ℓ M o ), and 0.0505 − 0.0928 pN/nm at 4 µm (8 µm or 3.25ℓ M o ).Since there are 6 titins acting in parallel for each half-sarcomere, we end up with the total stiffness for titin ranging between 0.0348 − 0.173 pN/nm at 2 µm and 0.303 − 0.557 pN/nm at 4 µm.When activated, the stiffness of our rabbit psoas linear-titin model (described in Sec.3.3, fitted in Appendix B.3, and with the parameters shown in Appendix H) doubles, which would increase titins stiffness to 0.0696 − 0.346 pN/nm at 2 µm and 0.606 − 1.11 pN/nm at 4 µm.

Figure 14 :
Figure14: The stiffness of a rabbit's actin-myosin load path with a single attached crossbridge (1 XB) exceeds the stiffness of its titin filament at lengths of 2 µm (1.62ℓ M o ) (compare AM:Low to TP:Low and TP:High).Only when titin is stretched to 4 µm (3.25ℓ M o ) does its stiffness (TP:High and TA:High) become comparable to the actin-myosin with a single attached crossbridge (AM:Low).At higher activations and modest lengths, the stiffness of the actin-myosin load path (AM: High) exceeds the stiffness of titin (TP: Low and TA:Low) by between two and three orders of magnitude.At higher activations and longer lengths, the stiffness of the actin-myosin load path (AM: High) exceeds the stiffness of titin by roughly an order of magnitude (TP:High and TA:High).
we can solve for k M in terms of k M AT by solving Eqn.36 for k M and substituting the values of Eqns.39, and 41.In this 6 case, the values of k M (4.46 N/mm) and k M AT (4.46 N/mm) B.3 Fitting the force-length curves of titin's seg-19 ments labelled the PEVK region of titin with antibodies allowing them to measure the distance between the Z-line and the proximal Ig/PEVK boundary ( Z ℓ IgP/PEVK ), and the distance between the Z-line and the PEVK/distal Ig boundary ( Z ℓ PEVK/IgD ), while the passive sarcomere was stretched from 2.35 − 4.46µm.By fitting functions to Trombitás et al.'s[26] data we can predict the length of any of titin's segments under the following assumptions: the T12 segment is rigid (Fig.1A), the distal16 Figure8of Prado et al.[58]  shows titin's contribution ranging from values ranging from (24%-57%) which means that the ECM's contribution ranges from (43%-76%) as a function of sarcomere length.Next, we extract the 17 coefficients for linear functions that evaluate the lengths of 18 lIgP ( lM ) = A IgP lM + b IgP , (53) lPEVK ( lM ) = A PEVK lM + b PEVK , and (54) lIgD ( lM ) = A IgD lM + b IgD (55)

18 Figure 12 22 Figure 3
Figure 12 data.The damping coefficients of the elastic and19

3 8 Y
can be transformed into the frequency-domain[78] by tak-4 ing the Fourier transform F(•) of the time-domain signal, 5 which produces a complex (with real and imaginary parts) 6 signal.Since convolution in the time-domain corresponds 7 to multiplication in the frequency-domain, we have (s) = H(s)X(s).

Figure 16 :
Figure 16: When coupled with a rigid-tendon, the VEXAT model's VAF (A.), gain response (B.), and phase response (C.) more closely follows the data of Kirsch et al. (Figure3)[5] than when an elastic-tendon is used.This improvement in accuracy is also observed at the 90 Hz perturbation (D., E., and F.), though the phase response of the model departs from Kirsch et al.'s data[5] for frequencies lower than 30 Hz. Parts of the Hill model's response to the 15 Hz perturbation are better with a rigid-tendon, with a higher VAF (G.), a lower RMSE gain-response (H.).but have a poor phase-response (I.).In response to the higher frequency perturbations, the Hill model's response is poor with an elastic (see Fig.6) or rigid-tendon.The VAF in response to the 90 Hz perturbation remains low (J.), and neither the gain (K.) nor the phase response of the Hill model (L.) follow the data of Kirsch et al.[5].The rigid-tendon Hill model's nonlinearity was so strong that the lowest frequency analyzed had to be raised from 4 Hz to 21 Hz to meet the criteria that (C xy ) 2 ≥ 0.67.

Figure 17 :
Figure 17: Simulation results of the 3 mm/s (A.) active lengthening experiment of Herzog and Leonard [7] (B.).As with the 9 mm/s trial, the Hill model's force response drops during the ramp due to a small region of negative stiffness introduced by the descending limb of the forcelength curve (C.), and a reduction in damping (D.) due to the flattening of the force-velocity curve.Note: neither model was fitted to this trial.

Figure 18 :
Figure 18: Simulation results of the 27 mm/s (A.) active lengthening experiment of Herzog and Leonard [7] (B.).As with the prior simulations the Hill model exhibits a small region of negative stiffness introduced by the descending limb of the force-length curve (C.) and a drop in damping (D.).Note: neither model was fitted to this trial.

6 segment bonds with actin through an activation-dependent 7 damper. The VEXAT model's distal segment of titin, ℓ 2 , 8 can contribute to RFE when the titin-actin bond is formed
11both our model and a Hill model remain ∆f B below the 12 experimental data of Herzog and Leonard (Fig.9B) after 13 lengthening and activation have ceased.
9an activation-dependent damper, when activation goes to 10 zero our titin model becomes unbound from actin.As such, 14 3.3 Active lengthening beyond actin-myosin over-15 lap 16 [8]nard 19et al.'s[8]active-lengthening experiment provides some 20 insight into this force distribution problem because they 21 recorded active forces both within and far beyond actin- 22myosin overlap.Leonard et al.'s [8]data shows that active 23 force continues to develop linearly during lengthening, be-24 Exp.
[7]have had to rely on the literature mentioned in Tablesimultaneouslyfit the passive force-length curve recorded 14 by Herzog and Leonard[7], using a mixture of tension from titin and the ECM that is consistent with Prado et al.'s data [58], all while maintaining the geometric relationship between f IgP and f PEVK as measured byTrombitás et al.
did 8 not report the architectural properties of the experimental 9 specimens, we assume that the architectural properties of 10 the cat used in Kirsch et al.'s experiments are similar to 11 the properties listed in Table 1.We evaluate the stiffness 12 of the tendon model by stretching it until it develops the 13 nominal tension of Kirsch et al.'s Figure 3 data (5N), and 14 then compute its derivative which amounts to k T =16.9 15 N/mm.Finally, using Eqn.33 we can solve for a value of 16 k M AT =2.90 N/mm.Since the inverse of damping adds for 17 damping elements in series

20
(50]e functions can be scaled to fit a titin filament of 21 a differing geometry.Many of the experiments simulated 22 in this work used cat soleus.Although the lengths of titin 23 filament segments in cat soleus have not been measured,24we assume that it is a scaled version of a human soleus titin 25 filament (68 proximal Ig domains, 2174 PEVK residues, 26 and 22 distal Ig domains[26]) since both muscles con-27 tain predominately slow-twich fibers: slow twitch fibers 28 tend to express longer, more compliant titin filaments[58].29Since the optimal sarcomere length in cat skeletal muscle 30 is shorter than in human skeletal muscle (2.43 µm vs. 2.73 31 µm,[75]) the coefficients for Eqns.53-55differslightly32(see the feline soleus column in Table2).In addition, by 33 assuming that the titin filament of cat skeletal muscle is a 34 scaled version of the titin filament found in human skeletal 35 muscle, we have implicitly assumed that the cat's skeletal 36 muscle titin filament has 60.5 proximal Ig domains, 1934.7 37 PEVK residues, and 19.6 distal Ig domains.Although a 38 fraction of a domain does not make physical sense, we have not rounded to the nearest domain and residue to preserve the sarcomere length-based scaling.In contrast, the rabbit psoas fibril used in the simulation of Leonard et al.[8]has a known titin geometry(50

Table 2 :
[53][54][55]ents of the normalized lengths of lIgP ( lM ), lPEVK ( lM ), and lIgD ( lM ) from Eqns.[53][54][55] [89]e the sarcomeres were incapable of developing active 10 or passive tension when the experiment was repeated after 11 the titin filaments were chemically cut.It is worth noting 12 that the forces measured by Leonard et al.[8]contain none 13 of the complex saw-tooth pattern indicative of unfolding 14 events even though 72 of these events would occur as each 15 proximal and distal Ig domain fully unfolded and reached 16 its maximal length17.Although we cannot be sure how 17 many unfolding events occurred during Leonard et al.'s ex-18 periments[8], due to sarcomere non-homogeneity[89], it 19 is likely that many Ig unfolding events occurred since the 20 19Using Eqns.53-55 and the appropriate coefficients in Table202 we can evaluate the normalized length developed by the 21 ℓ 1 segment CE length of lM * .The f 1 ( l 1 ) curve is formed by 24 shifting and scaling the (1−P)f PE curve so that it develops a normalized tension of (1 − P) and a stiffness of k 1 toe at a length of l1 toe .Similarly, the f 2 ( l 2 ) curve is made by shifting and scaling the (1 − P)f PE curve to develop a toe .
process we have used to fit the ECM and titin's seg- ing injury or during Leonard et al.'s experiment [8], 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 f 1 ( l 1 ) and f 2 ( l 2 ): a linear extrapolation, and the WLC model.Since the fitted f PE [5] 28[8].With the passive curves established, we can return to the 31 problem of identifying the normalized maximum stiffness 32 k X o and damping β X o of the lumped XE element.Just prior 33 to discussing titin, we had evaluated the impedance of the 34 cat soleus CE in Kirsch et al.'s[5]Figure12to be k M =2.90 35 N/mm and β M =0.020 Ns/mm at a nominal active tension 36 of 5N.The normalized stiffness k M can be found by taking 37 the partial derivative of Eqn. 15 with respect to lM k

Table 3 :
Normalized titin and crossbridge parameters fit to data from the literature.
T ).If we approximate the muscle's response as a linear 39( f we are interested in solving for H(s).While it the result will poorly estimate H(s) because Y (s) is only 11 approximated by H(s) X(s): Y (s) may contain nonlin-12 earities, non-stationary signals, and noise that cannot be 13 described by H(s) X(s).When the order of x(t) and y(t) are reversed 21 in Eqn.90 the result is G yx , while G xx and G yy are pro-22 duced by taking the Fourier transform of x(t) ⋆ x(t) and 23 y(t) ⋆ y(t) respectively.

Table 4 :
[5]h K segments allows us to evaluate lower frequency resolution than Eqn.89, but 32 an improved accuracy in H(s).Now we can evaluate the 33 gain of H(s) as Mean normalized stiffness coefficients (A.), mean normalized damping coefficients (B.), VAF (C.), and the bandwidth (D.) of linearity (coherence squared > 0.67) for models with elastic tendons.Here the proposed model has been fitted to Figure12of Kirsch et al.[5], while the experimental data from Kirsch et al.[5]comes from Figures 9 and 10.Experimental data from Figure 12 from Kirsch et al. has not been included in this table because it would only contribute 1 entry and would overwrite values from Figures9 and 10.The impedance experiments at each combination of perturbation amplitude and frequency have been evaluated at 10 different nominal forces linearly spaced between 2.5N and 11.5N.The results presented in the table are the mean values of these ten simulations.The VAF is evaluated between the model and the spring-damper of best fit, rather than to the response of biological muscle (which was not published by Kirsch et al.[5]).Finally, model values for the VAF (C.) and the bandwidth of linearity (D.) that are worse than those published by Kirsch et al.[5]appear in bold font.0.0118 0.0049 0.0038 0.0055 0.0041 0.0037 0.0198 0.0105 0.0024 0.8 mm 0.0118 0.0049 0.0038 0.0055 0.0042 0.0037 0.0165 0.0098 0.0027 1.6 mm 0.0118 0.0049 0.0038 0.0058 0.0043 0.0038 0.0122 0.0077 0.0023 34|H(s)| = (R(H(s)) 2 + I(H(s)) 2 ) (92)E Simulation summary data of Kirsch et al.24

Table 6 :
[8] VEXAT and Hill model's fitted rabbit psoas fibril MTU parameters.As in Table6, parameters shared by the VEXAT and Hill model are highlighted in grey.Short forms are used to save space: length 'len', velocity 'vel', acceleration 'acc', half 'h', activation 'act', segment 'seg', threshold 'thr', and stiffness 'stiff'.The letter preceding a reference indicates the experimental animal:'C' for cat, 'H' for human, while nothing at all is rabbit skeletal muscle.Letters following a reference indicate how the data was used to evaluate the parameter: 'A' for arbitrary for simulating Leonard et al.[8], 'n/a' for a parameter that is not applicable to a fibril model, '-' value taken from the cat soleus MTU, 'C' calculated, 'F' fit, 'E' estimated, 'S' scaled, and 'D' for default if a default value from another model was used.Only parameters that do not affect the outcome of our simulation of Leonard et al.[8]are marked 'A'.Clearly the parameters that appear in this Table do not represent a generic rabbit psoas fibril model, but instead a rabbit psoas fibril model that is sufficient to simulate the experiment of Leonard et al.[8].Finally, values for N IgP , N PEVK , and N IgD were obtained by taking a 70% and 30% average of the values for 3300 kD and 3400 kD titin to match the composition of rabbit psoas titin as closely as possible.