Understanding muscle function during perturbed in vivo locomotion using a muscle avatar approach

Decades of research on isolated muscles, single fibers and motor proteins has led to a widely accepted understanding of muscle mechanics based on the sliding filament–swinging cross-bridge theory (Huxley and Hanson, 1954; Huxley and Niedergerke, 1954; Rayment et al., 1993). Early studies emphasized the isometric force–length relationship (Gordon et al., 1966) and the isotonic force–velocity relationship (Hill, 1938), which was shown to be consistent with cross-bridge models of muscle contraction (Huxley, 1957, 1973). Compared with these traditional ex vivo methods for investigating muscle mechanics (i.e. isometric force–length and isotonic force–velocity relationships; Caiozzo, 2002), the work loop technique (Josephson, 1985) revolutionized the study of ex vivo muscle function (Ahn, 2012) by stimulating muscles during sinusoidal length oscillations at a frequency, phase and duty cycle similar to conditions that muscles experience during in vivo locomotion (Josephson, 1985, 1999). These studies enabled estimation of the work and power that muscles produce under in vivo conditions, compared with their maximum work and power output (Askew and Marsh, 1997; Askew et al., 1997; Full et al., 1998; Sponberg and Full, 2008; Sponberg et al., 2011; Tytell et al., 2018).

Work loop experiments also led to the insight that muscle function depends on the phase of activation relative to oscillations in length (Dickinson et al., 2000; Ahn et al., 2003). Muscles function like motors that perform positive work when stimulated at the onset of shortening, but like brakes that absorb work when stimulated during lengthening (Hessel and Nishikawa, 2017), or as struts and springs, producing little net work when stimulated during intermediate phases of the strain cycle (Ahn et al., 2003). These studies suggested that the role a muscle performs is determined largely by the timing of neural activation (Ahn et al., 2003; Ahn, 2012; Sponberg and Daniel, 2012).

Because of the technical difficulty of the measurements, relatively few studies have investigated muscle mechanics in vivo. In vivo animal studies have used sonomicrometry crystals to measure fascicle strain, electromyography (EMG) to measure activation, and tendon buckles to measure muscle force (Biewener et al., 1998; Daley et al., 2009; Daley and Biewener, 2011; Eng et al., 2019; Lee et al., 2013; Gordon et al., 2020; Roberts et al., 1997; Wakeling et al., 2021). These measurements are more difficult to make in humans, so muscle forces are typically estimated using inverse dynamics (Sylvester et al., 2021) rather than measured directly. However, a few recent studies have used ultrasound in humans (Dick et al., 2017) to measure not only muscle strain but also tendon strain, which should be proportional to muscle force.

Hill-type muscle models based on force–length and force–velocity relationships (Zajac, 1989; Seth et al., 2011) have been used to predict ex vivo (Stevens, 1993; James et al., 1995, 1996; Askew and Marsh, 1997; Perreault et al., 2003) and in vivo muscle forces (Lee et al., 2013; Dick et al., 2017; Wakeling et al., 2021). Yet, the results of these studies show that the force predictions of Hill-type models are relatively inaccurate (R²<0.6), especially at higher frequencies (James et al., 1996; Wakeling et al., 2021). The inability of Hill-type models based on the isometric force–length and isotonic force–velocity relationships to accurately predict muscle force in ex vivo experiments and during in vivo locomotion suggests that our understanding of muscle mechanics is incomplete (Umberger and Rubenson, 2011), and that some concept(s) or element(s) may be missing from the model.

One major limitation of both sinusoidal work loops and the force–velocity relationship is that they fail to include variability in muscle strain and loading rates that is typically observed during in vivo locomotion (Daley and Biewener, 2011). No current theory predicts muscle force during active shortening under time-varying loads and this is an important knowledge gap in muscle mechanics (Sponberg et al., 2023). Several new methods have been developed recently to overcome this limitation. For example, a novel bio-robotics approach has been used to apply realistic hydrodynamic loads to muscles using real-time feedback control, leaving the strain trajectory as the unconstrained and measured output variable (Richards, 2011; Richards and Clemente, 2012; Clemente and Richards, 2012; Robertson and Sawicki, 2015; Robertson et al., 2017). A similar approach was used by Olberding et al. (2019), who developed a computer model of limb extension during jumping to provide more realistic time-varying loads in ex vivo work loop experiments. However, these approaches are limited to relatively simple variations in loading that are amenable to either experimental manipulation or computational modeling. It remains to be discovered how naturally occurring variability in strain rates affects muscle force and work.

Technical advances have made it possible to measure in vivo strain and activation of individual muscles across a wide range of species and movement conditions, although direct measurements of force remain challenging for some muscles and species. To improve our understanding of variability in muscle function during in vivo conditions including perturbed locomotion, we developed a new method to bridge the gap between traditional ex vivo sinusoidal work loop experiments and in vivo locomotion. This method uses a well characterized muscle (extensor digitorum longus, EDL) from a laboratory rodent (mouse) as an ‘avatar’ to represent muscles that are difficult or impossible to study in ex vivo experiments. The method is otherwise similar to previous approaches using ex vivo work loop experiments to emulate in vivo strain and activation (e.g. Askew and Marsh, 1998, 2001; Sponberg et al., 2011). Development and validation of this avatar method enables testing the role that dynamic changes in strain and loading play during in vivo movements across a wide range of animals and muscles. We chose the term ‘avatar’ to describe use of the mouse EDL muscle to represent a different muscle in the ex vivo experimental environment, in which length and activation can be precisely controlled and the resulting forces measured.

The goal of the present study was to develop and validate the avatar method, and use it to investigate the contributions of naturally occurring small variations in muscle length and activation patterns to force production, measured in guinea fowl during a single trial of treadmill running with obstacle perturbations (Daley and Biewener, 2011). These obstacle perturbations produced stride-to-stride variation in muscle length, activation, force and work per cycle. To investigate mechanisms for stride to stride variation, we used mouse EDL muscles in ex vivo work loop experiments as an avatar for in vivo guinea fowl LG muscles. Our study used four strain trajectories from a single trial of perturbed running on a treadmill over obstacles at a speed of 2 m s⁻¹. We included strides up onto and down from obstacles attached to the treadmill, as well as unperturbed strides between obstacles. For comparison, we also included a purely sinusoidal strain trajectory at the same amplitude and frequency. This experimental design using ex vivo work loops allowed us to investigate how muscle intrinsic properties and stimulation interact to produce the variability in muscle force observed during perturbed in vivo locomotion. We also evaluated the ability of traditional sinusoidal work loops at the same frequency to predict muscle force. We used a combination of five strain trajectories and three stimulation patterns to examine the effects of variation in strain, stimulation and their interaction on muscle force and work in ex vivo work loops.

Based on Daley and Biewener (2011), we expected that muscle force production during active shortening would be associated not with length or velocity per se, but with abrupt changes in strain rate (i.e. strain transients) similar to those that result in vivo from variable loading associated with foot contact in different strides. Therefore, we predicted that in vivo strain trajectories from guinea fowl LG imposed on ex vivo mouse EDL muscles in avatar work loops would better explain stride-to-stride variations in muscle force than sinusoidal work loops that lack strain transients at the same amplitude and frequency. In purely isometric contractions, muscle force is related to activation, but during unexpected perturbations, muscles respond to changes in load faster than reflexes can modulate activation (Nichols and Houk, 1976). Therefore, we investigated the relative contributions of strain trajectories, stimulation patterns and their interaction on the timing and magnitude of peak muscle force and work per cycle in avatar work loop experiments.

EDL muscles (see fig. 1 of Chleboun et al., 1997) from 60–130 day old male wild-type mice (Mus musculus, n=6 muscles, each from a different mouse) of the strain B6C3Fe a/a-Ttn^mdm/J were used in this study. Breeder mice were obtained from the Jackson Laboratory (Bar Harbor, ME, USA) and a colony was established at Northern Arizona University (NAU). Mice were fed ad libitum and euthanized just prior to muscle extraction. The Institutional Animal Care and Use Committee at NAU approved the use of animals and experimental protocol.

EDL muscles were surgically exposed and 4-0 silk sutures were tied in square knots at the distal and proximal muscle–tendon junctions. Tendons were cut outside the suture knots and muscles were removed for ex vivo experiments. Extracted EDL muscles were attached to a dual-mode muscle lever system (Aurora Scientific, Inc., Series 300B, Aurora, ON, Canada) via tightened slip knots in the suture at the ends of the muscle. During experiments, the muscles were bathed in a 21°C, 7.4 pH Krebs–Henseleit solution containing (in mmol l^–1): NaCl 118, KCl 4.75, MgSO₄ 1.18, NaHCO₃ 24.8, KH₂PO₄ 1.18, CaCl₂ 2.54 and glucose 10.0. The bath was aerated with a 95% O₂/5% CO₂ gas mixture. The muscle was suspended between two platinum electrodes which delivered 1 ms square-wave pulses at tetanic supramaximal stimulation (80 V, 200 Hz) while finding optimal length (L₀). Submaximal stimulation (45 V, 90 Hz) was used during the experimental protocol to more closely emulate in vivo activation, and to decrease fatigue during experimental trials (James et al., 1995).

After finding L₀, a series of 80 V conditioning twitches was delivered to the muscle until twitch force reached a steady state (Hakim et al., 2013). Isometric force after a 500 ms tetanus was measured at L₀ using submaximal stimulation before and after the experimental protocol to measure fatigue. Muscles that lost more than 10% force (n=1) were not included in the analysis. After all experimental trials were completed, the muscle was removed from the rig, the sutures were removed, and the muscle was dabbed dry once and weighed. Physiological cross-sectional area (PCSA) was calculated using the standard formula: muscle mass/(L₀×1.06 g cm⁻³).

In a previous study (Daley and Biewener, 2011), strain trajectories, activation and force were measured in vivo from the lateral gastrocnemius (LG) of guinea fowl during running with obstacle perturbations on a treadmill, using sonomicrometry, EMG and tendon buckles (see fig. 1 of Daley and Biewener, 2003). Obstacle perturbations induced high stride-to-stride variability in muscle strain, activation and force (Fig. 1A). The guinea fowl LG muscle produced up to 40 times more work per cycle when the bird stepped up onto the obstacle than when it stepped down from the obstacle (Fig. 1A, compare red and gold lines). Additionally, the delay between EMG onset and force onset varied from as little as 20 ms for steps up onto the obstacle (Fig. 1A, red line) to more than 150 ms for steps down from the obstacle (Fig. 1A, gold line).

In vivo strain trajectories obtained from a single bird (Bl3; see Dryad dataset https://doi.org/10.7280/D11H49) running over obstacles at a speed of 2.0 m s⁻¹ (Fig. 1A) were used to control the length of the mouse EDL muscles during ex vivo work loop experiments (Fig. 1B). To investigate the function of the guinea fowl LG during different types of strides, we used LG strain trajectories from four strides of a guinea fowl's instrumented leg in ex vivo EDL work loop experiments: one stride in which the leg stepped down from the obstacle onto the treadmill (Fig. 1A; ‘Down’, gold line); two in which the instrumented leg stepped up from the treadmill onto the obstacle (Fig. 1A; ‘Up1’, purple line, and ‘Up2’, red line); and one in which no obstacle was present (Fig. 1A, ‘No obstacle’, blue line). We included two different Up trajectories because Up1 had a small stretch in the strain trajectory, whereas Up2 did not (Fig. 1A). Each of the strain trajectories was used in ex vivo EDL work loops at 2 Hz, although the in vivo guinea fowl strain trajectories varied in frequency from 2.85 to 3.34 Hz. The frequency was kept constant in the ex vivo experiments to isolate the effect of strain trajectory. Because in vivo data were sampled at a higher rate (10,000 Hz) than ex vivo data (4000 Hz), the in vivo data were down-sampled using MATLAB so that EDL strain, stimulation and force (Fig. 1B) were represented by 2000 data points per cycle. For each muscle, we also included passive and active sinusoidal work loops at the same strain amplitude and frequency. Two work loop cycles were performed for each strain trajectory. In preliminary studies (Rice, 2020), we tested mouse EDL muscles at different starting lengths (i.e. −10% L₀, −5% L₀, L₀, L₀+5%) and stretch amplitudes (10–20% L₀). Preliminary experiments using the Up1 strain trajectory suggested that an initial length of L₀−5% and a final length of L₀+5% best emulated the rate of increase in tension observed during passive stretch of the in vivo guinea fowl LG from 83.4% to 122.9% L₀ (see Table 1 and Fig. 1A).

Stimulation patterns for the ex vivo EDL work loops were based on measured in vivo guinea fowl EMG, which typically starts at the longest muscle length at the onset of leg retraction (Daley and Biewener, 2011), but starts ∼30 ms later in some birds (e.g. Ye3 in Dryad dataset https://doi.org/10.7280/D11H49). To measure effects of variation in stimulation patterns on work loop variables, all five strain trajectories were tested using three different stimulation patterns (EMG-based, Late and Long) at the same submaximal stimulation level (45 V, 90 Hz). On average, submaximal stimulation produced 91% of the maximum isometric force observed during supramaximal stimulation. At stimulation levels lower than 45 V and 90 Hz, the force in sinusoidal work loops was unfused. The onset and duration of the EMG-based stimulation pattern (Fig. 2A) used in the ex vivo EDL work loops was based on observed EMG data from the guinea fowl LG for the Up1 step. We accounted for the differing delay between activation and force onset, measured at ∼25 ms in vivo (Daley and Biewener, 2011) versus only 4–5 ms in ex vivo experiments, by stimulating EDL muscles 20 ms later than observed in vivo. Thus, EMG-based stimulation began 20 ms after the longest muscle length and continued for a duration of 115 ms as observed in vivo during the Up1 step (Daley and Biewener, 2011).

Late stimulation (Fig. 2B) was based on observed variation in EMG onset among birds (Daley and Biewener, 2011) and started 32.5 ms later than for EMG-based stimulation, for the same duration of 115 ms, occurring at the onset of increased muscle shortening velocity. Long stimulation (Fig. 2C) spanned the other two, as the muscle was stimulated at the same time as EMG-based stimulation, but was deactivated at the same time as Late stimulation, for an average duration of 148.4 ms (range 122 ms for Up1 to 156.25 ms for Down). Because the timing of stimulation onset was controlled relative to muscle length rather than cycle phase, the average duration varied among the strain trajectories. The order in which strain trajectories and stimulation patterns were performed was randomized for each muscle. In total, 15 work loops were performed on each muscle, including all combinations of five strain trajectories and three stimulation patterns.

The coefficient of determination (R²) was used to quantify the proportion of variance of in vivo guinea fowl LG force that was explained by the ex vivo mouse EDL force using the same strain trajectories. All 2000 data points per cycle were used to calculate R² values. The data were highly consistent within and between muscles. Given the observed standard deviations, the sample size of n=6 mice per group was sufficient to detect large effects (20% difference in group means) with power >0.8. The untransformed data for all five dependent variables (work per cycle; peak stress; muscle strain at peak stress; muscle strain at active stress onset; and time from muscle stimulation to peak stress) had non-homogeneous variances (Levene's test, P<0.05), and only peak stress was normally distributed (Shapiro–Wilk test, P<0.05). No transformations that we tried, including Box–Cox, succeeded in normalizing the data or equalizing the variances. Because the size of all model effects (stimulation pattern, strain trajectory and stimulation pattern×strain trajectory) was large for all dependent variables (P<0.008–0.0001), and because we were interested in the relative size of the main effects and their interaction, we conducted a three-way factorial analysis of variance despite violation of assumptions. ANOVA is generally robust for analysis of non-normal data even when sample sizes are small (Blanca et al., 2017), and also robust for analysis of heteroscedastic data when the sample sizes are equal (Glass et al., 1972). We also checked for non-independence of the five dependent variables using principal component analysis. When dependent variables are highly intercorrelated, the first principal component typically explains a high proportion (>75%) of the variance. For the work loop variables considered here, principal components 1–5 explained 39.8%, 25.7%, 17.3%, 14.1% and 3.1% of the total variance, respectively, demonstrating relatively low levels of intercorrelation.

Three-way mixed-model ANOVA was performed using JMP (version 14, SAS Institute Inc., Cary, NC, USA) to determine how stimulation (EMG-based, Late, Long; d.f.=2), strain trajectory (Up1, Up2, Down, No obstacle and Sinusoidal; d.f.=4) and their interaction affected work loop variables (peak stress, work per cycle, muscle strain at peak stress, muscle strain at force onset and time from stimulation onset to peak stress). The main effects were strain trajectory (fixed), stimulation pattern (fixed) and animal (random; d.f.=5). The interactions included animal×strain trajectory (random), animal×stimulation pattern (random), strain trajectory × stimulation (fixed) and animal×strain trajectory×stimulation pattern (random). Random effects, while sometimes significant, were not reported in ANOVA tables because they were not of interest in this study. To measure the relative size of the fixed effects (strain trajectory, stimulation pattern and strain trajectory×stimulation pattern), we analyzed FDR (false discovery rate) logWorth values. Tukey's honestly significant difference (HSD) test was used to evaluate post hoc differences among means.

As all ANOVA model effects were highly significant (P<0.008), we used discriminant function analysis to determine the relative importance of variation in strain trajectory and stimulation pattern in determining work loop variables. Linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA) were performed using JMP to determine whether trials could be discriminated based on the strain trajectories or stimulation patterns from which they were generated using the same dependent variables included in the ANOVA. These analyses were used to determine whether strain trajectory or stimulation pattern causes more distinct grouping of the data. A 50:50 test-to-train ratio was used in which the second of two work loops for each combination of strain trajectory and stimulation for each muscle was used to train the model and the first was used to test the model. Each strain trajectory group included a sample size of 18 work loops and each stimulation pattern group included a sample size of 30 work loops, representing trials with five strain trajectories and three stimulation patterns performed on six muscles for a total of 90 test trials for each discriminant analysis.

The mouse EDL and guinea fowl LG are both unipennate muscles (Table 1). Total muscle length was measured for the ex vivo mouse EDL muscle, but fascicle length was measured for the in vivo guinea fowl LG (Table 1). Because the pennation angle of the mouse EDL is small (12 deg), the difference between muscle length and fascicle length is also small (<3%). It is important to note, however, that the length and strain estimates from sonomicrometry may differ among locations within a given muscle (Higham and Biewener, 2011), and may not correspond to length measurements made using a servo-motor length sensor.

PCSA differed by a factor of 461 (from a mean of 0.0088 cm² in mouse EDL to 4.1 cm² in Bl3 LG), and the maximum observed forces differed by a factor of 177, from 0.25 N in mouse EDL (80.6% P₀) to 44.3 N (43.6% P₀) in the guinea fowl LG (Table 1). The observed maximum stress was 28.1 N cm⁻² for a strain amplitude of 10.3% L₀ in the mouse EDL versus 10.9 N cm⁻² for a strain amplitude of 39.5% median length in the guinea fowl LG. Therefore, the approximate active muscle ‘stiffness’ was ∼10 times higher for mouse EDL than for guinea fowl LG (Table 1).

Ex vivo strain trajectories were approximately sinusoidal (Fig. 1B), although they reached peak strain at different phases and were less symmetrical than the sinusoidal trajectory (Fig. 1B, black line). The average velocity was 0 mm s⁻¹ for all strain trajectories. A major difference between the sinusoidal and in vivo strain trajectories was that the in vivo strain trajectories contained abrupt changes in strain rate (i.e. strain transients) that were absent in the sinusoidal trajectory. During in vivo locomotion in the guinea fowl LG, these transients in strain rate correspond to foot contact (Daley and Biewener, 2011). The LG muscles typically reach their maximum shortening velocity immediately before foot contact and their maximum stretch velocity shortly thereafter. For the ex vivo experiments, the range of velocity was only 18.1 mm s⁻¹ for the sinusoidal trajectory, whereas it varied from 29.7 to 38.2 mm s⁻¹ among the in vivo strain trajectories.

When ex vivo mouse EDL work loops were scaled to the same approximate force and length as in vivo guinea fowl work loops (Fig. 3), the shapes of ex vivo EDL work loops were similar to the shapes of in vivo LG work loops with the same strain trajectory. R² values between in vivo guinea fowl LG force and ex vivo mouse EDL force ranged from 58% to 94%. EDL work loops with different strain trajectories showed distinctive features that were consistent across muscles; so much so, that work loops could be associated with a particular strain trajectory by looking at shape alone. The distinctive features of the work loops are likely due to strain transients that differed among the strain trajectories. The shapes of ex vivo mouse EDL work loops were more similar to in vivo guinea fowl LG work loops with the same strain trajectory than they were to ex vivo EDL work loops with different strain trajectories (Fig. 3).

In general, peak force, work per cycle, muscle strain at peak force and delay between activation/stimulation onset and peak force were similar for guinea fowl LG and mouse EDL, given the same strain trajectory (Fig. 3). The average variance in guinea fowl LG force explained (R²) by ex vivo mouse EDL force was 0.94 for the Up1 strain trajectory with Late stimulation, 0.76 for the Up2 strain trajectory with Late stimulation, 0.58 for the Down strain trajectory with Long stimulation, and 0.69 for the No obstacle strain trajectory with EMG-based stimulation. The variation in R² among the different strain trajectories is likely due to the fact that the stimulation patterns, amplitude, starting and ending strain values were optimized only for the Up1 strain trajectory, and those parameters were then applied to the other strain trajectories.

Despite similarities in overall work loop shape, there were differences between mouse EDL and guinea fowl LG (see Table 1). Although similar in length (fascicle L₀=14.6 mm for mouse EDL and 17.5 mm for guinea fowl LG), the mass, PCSA and maximum observed force of the guinea fowl LG were 627, 450, 461 and 177 times greater than those of the mouse EDL, respectively, and as noted previously, the mouse EDL was ∼10 times stiffer. Particularly for the Up1 trajectory, the mouse EDL failed to maintain peak force for as long as was observed in the guinea fowl LG (see Fig. 3A), suggesting that activation/deactivation kinetics are faster in the mouse EDL than in the guinea fowl LG.

In ex vivo work loop experiments, mouse EDL muscles exhibited considerable variation in force–length behavior among strain trajectories but with the same stimulation pattern (Fig. 4). The average R² was lower across strain trajectories for EDL muscles than within strain trajectories between mouse EDL and guinea fowl LG, keeping stimulation constant. For example, the average R² between mouse EDL Up1 Late and the other strain trajectories from the same EDL muscle was 0.29, compared with an R² of 0.94 between Up1 Late and the in vivo guinea fowl LG force from the same strain trajectory. The average R² values for mouse EDL Up1 Late compared with the other strain trajectories from the same muscle were as follows: Up2 Late R²=0.18, Down Long R²=0.33, No obstacle EMG-based R²=0.37. These comparisons suggest that ex vivo EDL stress was influenced more by strain trajectory than by muscle type (mouse EDL versus guinea fowl LG) or stimulation pattern (EMG-based, Long or Late).

Although the sinusoidal strain trajectory (Fig. 4, black line) looked somewhat similar to Up1 and Up2 strain trajectories and had the same strain amplitude (±5% L₀), frequency (2 Hz) and stimulation patterns, sinusoidal work loops were more rectangular in shape, as observed in previous studies (James et al., 1995, 1996). Regardless of stimulation pattern, sinusoidal strain trajectories lacked the distinguishing shape characteristics of ex vivo mouse EDL and in vivo guinea fowl LG work loops, and explained only a small proportion of variance in both guinea fowl LG forces for the same strain trajectories (average R²=0.045, n=4) and mouse EDL forces (average R²<0.109, n=6), which were nevertheless significantly greater than expected due to chance given the large number of data points (n=2000) on which the calculations were based.

Although stimulation patterns changed the shape of the EDL work loops, work loops with the same strain trajectory but different stimulation patterns were more similar to each other than work loops with the same stimulation pattern but different strain trajectories (Fig. 5). The effects of stimulation also differed among strain trajectories. All three stimulation patterns produced similar peak stress for the Up2 strain trajectory (Fig. 5B), but for the Down strain trajectory (Fig. 5C), EMG-based stimulation produced lower peak stress compared with the other stimulation patterns.

Strain trajectory, stimulation pattern and the strain×stimulation interaction had significant effects on peak stress (ANOVA, all P<0.0001; Table 2). No obstacle (mean±s.d. 25.7±0.60 N cm⁻²), Up1 (24.6±0.5 N cm⁻²) and Up2 (22.5±0.74 N cm⁻²) produced the highest peak stress, Down (21.3±1.2 N cm⁻²) produced significantly less stress than No obstacle and Up1, but was not statistically different from Up2, and the sinusoidal trajectory (12.8 ± 1.2 N cm⁻²) produced significantly lower peak stress than all other strain trajectories (Tukey HSD, all P<0.05; Fig. 6C). Late and Long stimulation produced significantly higher peak stress than EMG-based stimulation (Tukey HSD, all P<0.05). The interaction between strain trajectory and stimulation had a larger effect (FDR logWorth=24.4; Table 2) on peak stress than strain trajectory (FDR logWorth=7.7) or stimulation (FDR logWorth=7.2) alone. The main interaction effect was that all stimulation patterns produced nearly the same peak stress for Up1, Up2 and sinusoidal trajectories, but EMG-based stimulation produced significantly lower peak stress for the No obstacle and Down trajectories (Tukey HSD, all P<0.05; Fig. 7A).

Strain trajectory, stimulation and their interaction had significant effects on work per cycle (ANOVA, all P<0.0001; Table 2). Among strain trajectories, Up2 (115.9±7.1 mJ) and sinusoidal (89.6±10.5 mJ) produced significantly more work than Up1 (86.4±4.9 mJ), Down (86.6±7.3 mJ) and No obstacle (71.2±4.7 mJ; Tukey HSD, all P<0.05; Fig. 6A). Among activation patterns, Long stimulation produced significantly more work per cycle than the other stimulation patterns, and Late stimulation produced more work per cycle than EMG-based stimulation (Tukey HSD, all P<0.05). The main interaction effect was that EMG-based and Late stimulation produced nearly the same work per cycle for Up1, Up2 and sinusoidal strain trajectories, but Late stimulation produced more work per cycle than EMG-based stimulation for No obstacle and Down strain trajectories (Tukey HSD, all P<0.05; Fig. 6A). The interaction between strain trajectory and stimulation pattern had a larger effect on work per cycle (FDR logWorth=17.9) than either strain trajectory (9.4) or stimulation (4.6) alone (Table 2).

Strain trajectory, stimulation and the strain trajectory×stimulation interaction had significant effects on muscle strain at peak stress (ANOVA, all P<0.01; Table 2). Muscle strain at peak stress was significantly higher for Up1 (0.99±0.001% L₀) and Up2 (1.00±0.0004% L₀) compared with No obstacle (0.96±0.0002% L₀), Down (0.97±0.005% L₀) and Sine (0.98±0.002% L₀; Tukey HSD, all P<0.05; Fig. 6D). EMG-based stimulation produced significantly larger strains at peak stress than other stimulation patterns (Tukey HSD, all P<0.05). The main interaction effect was that all stimulation patterns produced nearly the same strain at peak stress for Up1, Up2 and No obstacle, but EMG-based stimulation produced significantly larger strains at peak stress for the Sinusoidal and Down trajectories (Tukey HSD, all P<0.05; Fig. 7A). Strain trajectory had a larger effect (FDR logWorth=7.7; Table 2) on strain at peak stress than stimulation (FDR logWorth=3.0), or strain trajectory×stimulation (FDR logWorth=2.1).

Strain trajectory, stimulation and their interaction had significant effects on muscle length at active stress onset (ANOVA, all P<0.005; Table 2). Sinusoidal (1.045±0.001% L₀), Up2 (1.044±0.001% L₀) and Down (1.044±0.004% L₀) trajectories developed stress at a significantly longer length than Up1 (1.037±0.002% L₀) and No obstacle (1.036±0.002% L₀; Tukey HSD, all P<0.05; Fig. 6A). Strain at stress onset was significantly lower with Late than with EMG-based or Long stimulation (Tukey HSD, all P<0.05). The main interaction effect was that EMG-based and Long stimulation had nearly the same strain for all trajectories except Down, where EMG-based stimulation produced significantly higher strain at stress onset than Long stimulation. Late stimulation produced significantly lower strain than the other stimulation patterns for Up1, Up2 and No obstacle trajectories, but was not statistically different from Long stimulation for the Down trajectory. All stimulation patterns produced nearly the same strain at stress onset for the Sinusoidal strain trajectory (Tukey HSD, all P<0.05; Fig. 6B). Stimulation had a larger effect (FDR logWorth=4.801; Table 2) on strain at stress onset than strain trajectory×stimulation (FDR logWorth=4.111) or strain trajectory (FDR logWorth=2.425).

Strain trajectory, stimulation and their interaction had significant effects on time from activation onset to peak stress (ANOVA, P=0.001; Table 2). No obstacle (136±2.4 ms), Sine (126±3.1 ms) and Down (131±6.5 ms) produced statistically similar times from stimulation to peak stress. Up1 (111±2.5 ms) produced the shortest time, but was not statistically different from Up2 (123±1.5 ms; Tukey HSD, all P<0.05; Fig. 6E). Long stimulation produced significantly longer times from stimulation to peak stress than other stimulation patterns (Tukey HSD, all P<0.05). The main interaction effect was that all stimulation patterns produced times that were not statistically different for the Up1 and Up2 strain trajectories. Long stimulation produced significantly longer times for No obstacle, Sinusoidal and Down trajectories. EMG-based and Late activation produced nearly the same times for all trajectories except Down, while Late stimulation produced significantly longer times (Tukey HSD, all P<0.05; Fig. 7A). Stimulation had a larger effect (FDR logWorth=5.006; Table 2) on time from stimulation to peak stress compared with strain trajectory×stimulation (FDR logWorth=4.080) and strain trajectory (FDR logWorth=3.837).

In summary, strain at stress onset and time from stimulation to peak stress were more affected by stimulation than by strain trajectory or the interaction between these variables. Strain at peak stress was most affected by strain trajectory, while work per cycle and peak stress were more affected by the strain trajectory×stimulation interaction.

Because strain trajectory, stimulation pattern and their interaction had highly significant effects (P<0.008) on all work loop variables, we used discriminant function analysis to quantify how the strain trajectories (Fig. 7A) and stimulation patterns (Fig. 7B) cluster in the space defined by the five work loop variables included in the ANOVA. QDA was more accurate than LDA in classifying the trials, likely because it allows for unequal covariances among groups. Using the work loop variables to categorize strain trajectories (Fig. 7A), QD1 and QD2 explained 98.2% of the total information (entropy R²), with length at peak stress and peak stress having the highest weightings. All strain trajectories were classified 100% correctly (n=18 for each strain trajectory).

Results of the QDA analysis to classify stimulation patterns (Fig. 7B) showed that QD1 and QD2 explained 74.5% of total information, with time from stimulation to peak stress and work per cycle having the highest weightings. Five Late stimulation trials were misclassified as Long, while nine Long stimulation trials were misclassified as Late (14 of 90 trials misclassified, 84.4% correct). The results from this dataset show that variation among trials in all work loop variables considered simultaneously was more closely related to strain trajectory than to stimulation pattern.

Compared with traditional ex vivo methods for investigating muscle mechanics (e.g. isometric force–length and isotonic force–velocity characteristics; Caiozzo, 2002), the work loop technique (Josephson, 1985) revolutionized the way ex vivo muscle function is investigated (Ahn, 2012), by stimulating muscles during sinusoidal oscillations in length at a frequency, phase and duty cycle similar to conditions that muscles experience during steady in vivo locomotion (Josephson, 1985, 1999). These studies enabled estimation of the work and power that muscles can produce under in vivo conditions, compared with their maximum power output (Askew and Marsh, 1997, 2001), and further showed that saw-tooth trajectories, with more variable strain rates than purely sinusoidal strain trajectories, produce more force, work and power (Marsh, 1999). Like the present study, many previous studies used ex vivo work loops with in vivo strain and activation; for example, sonomicrometry and EMG recordings from quail pectoralis during vertical take-off (Askew and Marsh, 2001). However, these studies were based on in vivo strain trajectories under steady conditions in which cycle-to-cycle variability in muscle strain is relatively small (see fig. 2B of Askew and Marsh, 2001).

The novelty of the avatar method is to use a widely studied and readily available muscle (i.e. mouse EDL) to represent a different muscle in ex vivo work loop experiments. The method can be used to investigate muscle mechanics under in vivo-like conditions for muscles and species that are otherwise unsuitable for ex vivo experimentation, including humans and other large animals. In the present study, we used the avatar method to investigate mechanisms that lead to the large variability in peak muscle force and work exhibited in vivo by the guinea fowl LG muscle during perturbed locomotion.

Perturbations provide insights into muscle intrinsic properties that are not evident under steady-state conditions (Daley and Biewener, 2011; Shorter and Rouse, 2018; Sponberg et al., 2023). This study used cyclical length changes in ex vivo experiments to understand how muscles function during perturbed in vivo locomotion. The results of this study demonstrate that: (1) ex vivo mouse EDL work loops explained a large proportion (58–94%) of the variance in the in vivo force of guinea fowl LG given similar strain and stimulation, despite many differences between muscles and conditions; (2) strain trajectory, stimulation pattern and strain trajectory×stimulation interaction had highly significant (P<0.008) effects on all work loop variables, and for peak stress and work per cycle the interaction effect was larger than the effects of either strain trajectory or stimulation pattern alone; and (3) in vivo strain trajectories produced consistent ex vivo work loops that could be identified with high accuracy (100%) based on work loop variables using discriminant function analysis. These results demonstrate that using perturbed strain trajectories derived from in vivo measurements significantly improves prediction of stride-to-stride variability in muscle force during in vivo perturbed locomotion, compared with sinusoidal strain trajectories that lack strain transients. Strain transients are associated with changes in muscle loading during in vivo locomotion (Daley and Biewener, 2011). Because it applies strictly to isotonic contractions, Hill's (1938) force–velocity relationship fails to predict muscle force under conditions of variable strain and load (Libby et al., 2020; Jeong and Nishikawa, 2023). To improve understanding of in vivo muscle function under dynamic loading conditions, more investigations that measure strain, activation and force of individual muscles during perturbed locomotion are needed, as well as new ex vivo techniques that measure the force response of muscles under varying loading and strain. These efforts should lead to muscle models that can more accurately predict in vivo muscle forces.

Ex vivo EDL experiments differed from in vivo guinea fowl experiments in several ways. During in vivo experiments, there was variation in EMG activity between steps, even between steps of the same type (Daley and Biewener, 2011). In contrast, our experiments used submaximal square-wave electrical stimulation (45 V, 90 Hz). Furthermore, the ex vivo and in vivo data were from different muscles (LG versus EDL) and animals (guinea fowl versus mouse). Despite these differences, in vivo strain trajectories from guinea fowl LG produced highly consistent work loops in mouse EDL with similar variability in peak stress, work per cycle, muscle strain at force onset and delay between stimulation onset and peak force. The average variance in guinea fowl LG force explained by EDL avatar force was 58–94%, with the highest R² value obtained for the trial with EMG-based stimulation. The EMG-based stimulation pattern produced significantly lower peak stress for the No obstacle and Down trajectories, which may help to explain why guinea fowl use the in vivo activation patterns that they do. No obstacle and Down strides involve more energy dissipation than the other step types, and the EMG-based stimulation pattern generated lower peak stress for these stride types, so that there was less energy to absorb.

The R² values from ex vivo EDL work loops are similar to or exceed Hill model predictions of in vivo muscle forces (Lee et al., 2013; Dick et al., 2017; Wakeling et al., 2021), suggesting that ex vivo work loop experiments using muscles from laboratory rodents might serve as an alternative ‘model’ for understanding in vivo muscle mechanics in general. Development of techniques for ex vivo muscle stimulation that enable more realistic time-varying amplitudes and frequencies that better emulate in vivo activation patterns is an important direction for future studies.

It is somewhat surprising that ex vivo work loops from different muscles and animals, with simplified electrical stimulation compared with in vivo neural activation, generate forces similar to the in vivo LG force when given similar strain trajectories. The overall similarities in force–length behavior between guinea fowl LG and mouse EDL demonstrate the importance of strain, and particularly strain transients, in muscle force production during in vivo locomotion. The similarities suggest that the intrinsic viscoelastic properties of muscles play an important role in force production by resisting abrupt changes in length, even with no stretching. The resulting forces are not represented in sinusoidal oscillations or predicted by the isotonic force–velocity relationship (Libby et al., 2020; Nishikawa et al., 2018; Jeong and Nishikawa, 2023).

The present study also included sinusoidal strain trajectories at the same amplitude (±5% L₀) and frequency as the in vivo strain trajectories. The sinusoidal strain trajectories produced rectangular-shaped work loops, as described in previous studies (James et al., 1995, 1996), that bore little resemblance (average R²=0.045) to the triangular and L-shaped work loops produced in vivo (see Fig. 3), despite overall similarity to the shapes of the in vivo strain trajectories (see Fig. 1B). Sinusoidal strain trajectories also produced significantly lower peak forces than the in vivo strain trajectories (see Fig. 6C). The main difference between the sinusoidal and in vivo strain trajectories is that the in vivo trajectories contain strain transients associated in vivo with foot contact that are strongly linked with muscle force onset (see Fig. 1B), whereas the sinusoidal strain trajectory at the same frequency does not.

The results suggest that using sinusoidal, saw-tooth or smoothed in vivo strain trajectories in ex vivo work loop experiments may result in underestimation of peak forces and may fail to accurately represent in vivo muscle mechanics. Many previous studies did not include strain transients in ex vivo work loop experiments to investigate in vivo muscle function (Askew and Marsh, 2001; Tytell et al., 2018), in part because they studied swimming and flying where changes in loading are less apparent than in terrestrial locomotion. However, strain transients are likely to occur during non-steady locomotion in any medium, and are expected to have large effects on in vivo muscle function.

Studies using the work loop technique have also demonstrated that, by changing the phase of stimulation relative to ongoing length changes, the function of a muscle can change from a motor to a strut, spring or brake (Dickinson et al., 2000; Ahn et al., 2003). These studies revealed the importance of the nervous system in controlling the mechanical output of muscles (Ahn, 2012). However, it is also necessary to consider the stabilizing function of muscle intrinsic properties during unexpected perturbations (i.e. ‘preflexes’; Loeb, 1995), which can modulate muscle mechanical output without requiring feedforward or feedback instructions from the nervous system. We found that the strain trajectory has a large influence on peak stress and work per cycle of ex vivo muscles. Although the guinea fowl LG is active during shortening and functions mostly like a motor during walking and running (Daley and Biewener, 2011), the ex vivo EDL data demonstrate that, given constant activation, strain trajectories strongly influence the peak force and work per cycle that a muscle can produce, increasing positive work during strides up onto an obstacle, and dissipating work during strides down from an obstacle, compared with steady locomotion.

The results of this study support the idea that muscle is an active material whose viscoelastic properties are tuned by activation (Kirsch et al., 1994; Nguyen and Venkadesan, 2020preprint; Nishikawa, 2020). The major difference between this view of muscle and the traditional view of muscle as a motor that produces force in response to activation is that information, not only about activation but also about strain transients, is necessary to predict the force of a tunable material. In ex vivo experiments, it is possible to examine the effects of different strain trajectories on force production while holding stimulation constant. Our results demonstrate that the force–length behavior of muscles in work loop experiments is more closely associated with strain trajectory than with stimulation pattern (see Fig. 7). Furthermore, work loop variables including peak force and work per cycle were more affected by the interaction between activation and strain trajectory than by either variable alone (see Fig. 6, Table 2). Only 2 of 5 work loop variables (length at stress onset and time from stimulation to peak force) were more strongly affected by stimulation than by strain trajectory or the interaction between these variables (Fig. 6, Table 2). The ex vivo results account for the observed decoupling in time between activation and force during perturbed locomotion of guinea fowl (Daley et al., 2009; Daley and Biewener, 2011), and support the idea that muscle is a tunable viscoelastic material for which activation permits but does not instruct force production. The onset and magnitude of muscle force production are determined not only by activation but also importantly by loads applied to the muscle during interaction with the environment (Nishikawa et al., 2007). Understanding the interaction between strain and activation is key to understanding in vivo muscle function.

In this study, the muscle mechanics data were highly consistent within and between muscles for all combinations of stimulation patterns and strain trajectories. Future experiments should include validating the avatar method using in situ guinea fowl LG muscles. A comparison of results from ex vivo mouse EDL muscles with those from in situ guinea fowl muscles would enable quantification of the relative contributions of muscle size, fiber type, activation and deactivation kinetics, etc., to muscle force generation. Future experiments should also explore development of techniques for ex vivo muscle stimulation that enable more realistic time-varying amplitudes and frequencies that better emulate in vivo activation patterns.

We thank Benjamin Thiesing for technical assistance. Anthony Hessel, Natalie Holt, Tinna Traustadóttir and Philip Service provided guidance and revisions on earlier versions of this manuscript.