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 the titin and actin in proportion to the activation of the muscle (C.), which reduces to negligible values in a purely passive muscle (D.). The ECM is the only entirely passive element (A. & D.).

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 a 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 [6]. 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 [60] (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. [24] (C.) observed in their experiments. Note that the curves that appear in subplots A. and B. come from scaling the generic titin model, described in Appendix B, to a cat soleus sarcomere while the curves that appear in subplot C come from scaling the generic titin model to a human soleus sarcomere to compare it to the data of Trombitás et al.’s [24].

Kirsch et al. [5] performed an in-situ experiment on a cat soleus. The soleus was activated to generate a constant submaximal tension at a fixed length. Next, the length of the musculotendon unit was forced to follow a stochastic waveform (see Fig. 6 for an example) while the force response of the muscle was recorded.

Different arrangements of springs and dampers will produce different force responses when subjected to the same length perturbation waveform in the time-domain (A.). In the time-domain it is challenging to identify the system given only the input length change and output force signals. If the system is assumed to be linear and time invariant then system identification methods [71] (see Appendix D) can be used to evaluate the gain (B.) and phase (C.) response using only the input length change and an output force profile. Given the gain (B.) and phase (C.) response it is often possible to identify the system due to the clearly identifiable responses of each system.

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 88-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 too steep (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 frequency coefficients of both the VEXAT model and the Hill model are scattered points (rather than continuous curves) because of the nonlinearities present in each model.

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 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 stiffness and damping coefficients of 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 31%-44% range shown in A. and B. to 79% shown in C. and D. Although the Hill model’s response more closely resembles a spring-damper when coupled with a rigid tendon, the spring and damper coefficients fail to fit Kirsch et al.’s data [5].

Kirsch et al. [5] noted that the VAF of the spring-damper model of best fit captured between 88-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 VAF of the VEXAT model is affected by up to 20% depending on the condition, with the lowest VAF occurring in response to the 15 Hz 1.6 mm perturbation (A.), with higher VAF responses occurring at 35 Hz (B.) and 90 Hz (C.). Where the tendon model affects the VAF of the VEXAT model slightly, the tendon model causes the VAF of the Hill model to vary by nearly 60% depending on the condition: at the 15 Hz (A.) and 35 Hz (B.) the Hill model’s response is dramatically affected by the tendon model, while at the 90 Hz (C.) the tendon model does not exert such a large influence.

Herzog and Leonard [6] actively lengthened (A.) a cat soleus 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. [7] 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 [8] 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 1.05 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 [9] 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.

Cat soleus MTU properties used in this work. Parameter Symbol Value Source

The normalized tendon force-length curve (A) has been been fit to match the cat soleus tendon stiffness measurements of Scott and Loeb [79]. 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).

The coefficients of the normalized lengths of , and from Eqns. 51-53 under passive lengthening. These coefficients have been extracted from data of Figure 5 of Trombitás et al. [24] using a least-squares fit. Since Figure 5 of Trombitás et al. [24] 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 sacromere (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 [60] differs drammatically from human soleus titin [24] 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.), damping coefficients (B.) and VAF (C.) for models with elastic tendons. Here the proposed model has been fitted to Figure 12 of Kirsch et al. [5]. The impedance experiments at each combination of perturbation amplitude and frequncy have been evaluated at 10 different nominal forces evenly spaced between 2.5N and 11.5N. The normalized results presented in the table are the mean values of these ten simulations. Finally, note that the VAF is evaluated between the model and the spring-damper of best fit to the response of the model, rather than to the response of biological muscle (which was not published by Kirsch et al. [5]).

Mean normalized stiffness coefficients (A.), damping coefficients (B.) and VAF (C.) for models with rigid tendons. All additional details are identical to those of Table 4 except the tendon of the model is rigid.

Simulation results of the 3 mm/s (A.) active lengthening experiment of Herzog and Leonard [6] (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 [6] (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.

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 (B.). but a similarly poor phase-response (C.). In response to the higher frequency perturbations, the Hill model’s response is poor with an elastic (see Fig. 7) or rigid tendon. The VAF in the time-domain remains low (J.), neither the gain (K.) nor the phase response of the Hill model (L.) follow the data of Kirsch et al. [5].