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

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

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

The 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 and . 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). Note that the rigid-tendon cat soleus parameters differ from the table below because tendon elasticity affects the fitting of , and ). Finally, the parameters related to the compressive element (F.), the XE (G.), and titin (H. and I.) can be used as initial values when simulating the MTU’s other mammals. By making use of these defaults the VEXAT model requires the same number of parameters as a Hill-type muscle model (A.—E.).

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

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.).

The 15 Hz perturbations show that the VEXAT model’s performance is mixed: in the time-domain (A.) the VAF is lower than the 78-99% analyzed by Kirsch et al. [5]; the gain response (B.) follows the profile in Figure 3 of Kirsch et al. [5], while the phase response differs (C.). The response of the VEXAT model to the 90 Hz perturbations is much better: a VAF of 91% is reached in the time-domain (D.), the gain response follows the response of the cat soleus analyzed by Kirsch et al. [5], while the phase-response follows biological muscle closely for frequencies higher than 30 Hz. Although the Hill’s time-domain response to the 15 Hz signal has a higher VAF than the VEXAT model (G.), the RMSE of the Hill model’s gain response (H.) and phase response (I.) shows it to be more in error than the VEXAT model. While the VEXAT model’s response improved in response to the 90 Hz perturbation, the Hill model’s response does not: the VAF of the time-domain response remains low (J.), neither the gain (K.) nor phase responses (L.) follow the data of Kirsch et al. [5]. Note that the Hill model’s 90 Hz response was so nonlinear that the lowest frequency analyzed had to be raised from 4 Hz to 7 Hz to satisfy the criteria that (Cxy)2 ≥ 0.67.

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 and have been calculated to take the tendon model into account (see Appendix B.4 for details). 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 increases the VAF from 30-51% to 86% but this improved accuracy is made using stiffness and damping that deviates from that of biological muscle [5].

Kirsch et al. [5] noted that the VAF of the spring-damper model of best fit captured between 78-99% across all experiments. We have repeated the perturbation experiments to evaluate the VAF across a range of conditions: two different tendon models, three perturbation bandwidths (15 Hz, 35 Hz, and 90 Hz), three perturbation magnitudes (0.4 mm, 0.8 mm and 1.6 mm), and ten nominal force values (spaced evenly between 2.5N and 11.5N). Each bar in the plot shows the mean VAF across all ten nominal force values, with the whiskers extending to the minimum and maximum value that occurred in each set. The mean VAF of the VEXAT model changes by up to 36% depending on the condition, with the lowest mean VAF occurring in response to the 1.6 mm 15 Hz perturbation with an elastic-tendon (A.), and the highest mean VAF occurring in response to the 90 Hz perturbations with the rigid-tendon (C.). In contrast, the mean VAF of the Hill model varies by up to 67% depending on the condition, with the lowest VAF occurring in the 15 Hz 0.4 mm trial with the elastic-tendon (A.), and the highest value VAF occurring in the 15 Hz 0.4 mm trial with the rigid-tendon (A.).

Herzog and Leonard [7] actively lengthened (A.) cat soleus muscles on the descending limb of the force-length curve (where 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.).

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 .

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

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 to the left the VEXAT model (magenta squares) can be made to exactly match the force-length characteristic of the Hill model.

To evaluate the stiffness of the actin-myosin load path, we first determine the average point of attachment. Since the actin filament length varies across species we label it LA. Across rabbits, cats and human skeletal muscle myosin geometry is consistent [75]: a half-myosin is 0.8µm in length with a 0.1µm bare patch in the middle. Thus at full overlap the average point of attachment is 0.45µm from the M-line, or LA−0.45µm from the Z-line at . 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 LA − 0.45µm length of actin, a spring representing the all attached crossbridges, and a spring representing a 0.45µm section of myosin.

The stiffness of a rabbit’s actin-myosin load path with a single attached crossbridge (1 XB) exceeds the stiff-ness of its titin filament at lengths of (com-pare AM:Low to TP:Low and TP:High). Only when titin is stretched to 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).

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).

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

The coefficients of the normalized lengths of , and from Eqns. 53-55 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.

Normalized titin and crossbridge parameters fit to data from the literature.

Mean normalized stiffness coefficients (A.), mean normalized damping coefficients (B.), VAF (C.), and the bandwidth (D.) of linearity (coherence squared > 0.67) for models with elastic-tendons. Here the proposed model has been fitted to Figure 12 of Kirsch et al. [5], while the experimental data from Kirsch et al. [5] comes from Figures 9 and 10. Experimental data from Figure 12 from Kirsch et al. has not been included in this table because it would only contribute 1 entry and would overwrite values from Figures 9 and 10. The impedance experiments at each combination of perturbation amplitude and frequency have been evaluated at 10 different nominal forces linearly spaced between 2.5N and 11.5N. The results presented in the table are the mean values of these ten simulations. The VAF is evaluated between the model and the spring-damper of best fit, rather than to the response of biological muscle (which was not published by Kirsch et al. [5]). Finally, model values for the VAF (C.) and the bandwidth of linearity (D.) that are worse than those published by Kirsch et al. [5] appear in bold font.

Mean normalized stiffness coefficients (A.), mean normalized damping coefficients (B.), VAF (C.), and the bandwidth (D.) of linearity (coherence squared > 0.67) for models with rigid tendons. All additional details are identical to those of Table except the tendon of the model is rigid.

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

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

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.