Adjustments of motor pattern for load compensation via modulated activations of muscle synergies during natural behaviors

Abstract

It has been suggested that the motor system may circumvent the difficulty of controlling many degrees of freedom in the musculoskeletal apparatus by generating motor outputs through a combination of discrete muscle synergies. How a discretely organized motor system compensates for diverse perturbations has remained elusive. Here, we investigate whether motor responses observed after an inertial-load perturbation can be generated by altering the recruitment of synergies normally used for constructing unperturbed movements. Electromyographic (EMG, 13 muscles) data were collected from the bullfrog hindlimb during natural behaviors before, during, and after the same limb was loaded by a weight attached to the calf. Kinematic analysis reveals the absence of aftereffect on load removal, suggesting that load-related EMG changes were results of immediate motor pattern adjustments. We then extracted synergies from EMGs using the nonnegative matrix factorization algorithm and developed a procedure for assessing the extent of synergy sharing across different loading conditions. Most synergies extracted were found to be activated in all loaded and unloaded conditions. However, for certain synergies, the amplitude, duration, and/or onset time of their activation bursts were up- or down-modulated during loading. Behavioral parameterizations reveal that load-related modulation of synergy activations depended on the behavioral variety (e.g., kick direction and amplitude) and the movement phase performed. Our results suggest that muscle synergies are robust across different dynamic conditions and immediate motor adjustments can be accomplished by modulating synergy activations. An appendix describes the novel procedure we developed, useful for discovering shared and specific features from multiple data sets.

FIG. 1.

Experimental procedure and the muscle synergy model used in this study. A: the experimental schedules of the 4 frogs (P, Q, R, and S) studied herein. Electromyographic data of 13 muscles during terrestrial behaviors (kicking, jumping, and stepping) were collected from 3 frogs (P, Q, and R) and those during aquatic behaviors (in-phase and out-of-phase swimming), from one frog (S). In selected trials, an inertial load was applied to the frog hindlimb by attaching a weight to the frog's calf. The black boxes in the figure indicate the within-day sessions in which each frog was exposed to the load. The number in each box refers to the weight of the load applied, expressed as a percentage of the hindlimb's weight. Abbreviations: load1, the first block of loaded trials; load2, the second block of loaded trials; load3, the third block of loaded trials. B: the muscle synergy model. In this example, the waveforms of 3 muscles (M1, M2, and M3) are explained as a linear combination of 2 synergies (synergy 1 and synergy 2), each represented as a static activation balance profile across the 3 muscles, and activated by a time-varying activation coefficient. The model can be formally stated by the equation shown in the figure, where D(t) is the electromyographic (EMG) data of all muscles at time t, w_i is the synergy vector for the ith synergy, c_i(t) is the coefficient for the ith synergy at time t, and M is the total number of synergies composing the data set. For more details, see methods, Data analysis: the muscle synergy model.

FIG. 2.

Characterizing varieties of kicks and jumps. A: as shown in the leftmost frog diagram, the hip angle is defined as the angle between the body axis and femur; the knee angle is defined as the angle between femur and tibiofibula. The rest of the frog diagrams in this panel illustrate the kick displacement vector, defined as the joint-angle difference between the initial ankle position (h1 and k1) and the kick's maximally extended position (h2 and k2). The hip and knee joint-angle displacements define the x- and y-axes of the kick behavioral space, respectively. B: the extents of distribution of kick varieties in the kick behavioral space before and after inertial loading. As can be seen here, the extents of distribution of unloaded (+) and loaded (○) kick varieties almost completely overlap each other. The 3 frog cartoons shown in the figure indicate where lateral kicks (top left), caudal–lateral kicks (top right), and caudal kicks (bottom) lie in the kick behavioral space. C: the distribution extents of the unloaded and loaded kicks were quantified by first defining a baseline set of kicks (for definition, see methods, Trajectory analysis) and then calculating the Euclidean distance between each kick and its closest baseline neighbor in the behavioral space. After binning the kicks by kick magnitude (defined as the magnitude of the kick displacement vector), we plotted the distances for the loaded (○) and unloaded (▪) kicks (mean ± SD). For all bins (except the first bin of frog Q, marked by *), the unloaded and loaded distances are not significantly different from each other (P > 0.05). This suggests that the differences between the loaded and baseline kicks are not more than the differences among the unloaded kicks. D: jump length is defined as the distance of the line connecting the initial and final midpoint positions of the frog. Jump direction is defined as the angle between the initial body axis and the line connecting the frog's initial and final midpoint positions. Positive angles are assigned to rightward jumps and negative angles, to leftward jumps. E: distribution extents of jump varieties before (+) and after (○) inertial loading. The jump behavioral space is represented by both jump length (y-axis, expressed as a multiple of the body length) and jump direction (x-axis). The 2 distributions almost completely overlap each other in extent.

FIG. 3.

The inertial load's effect on kick and jump trajectories. A and C: for each kick (A) and jump (C) episode, the average distance between the episode's joint-space trajectory (with space origin denoting the initial limb position) and those of the 5 baseline neighbors closest to that episode in the behavioral space was calculated. Trajectory distance was quantified by the figural distance measure introduced by Conditt et al. (1997). Plots of this distance against time (represented as chronologically ordered episode number) show that for both kicks (A) and jumps (C), there was no increase in figural distance in the beginning of the postload epochs, suggesting that EMG changes associated with loading represent immediate motor pattern adjustments rather than motor learning. B and D: in the loaded epochs (load1, load2, and load3), there were several episodes with notably higher figural distances (as shown in A and C). Shown here are positions of the loaded kicks (B) and loaded jumps (D) in their respective behavioral spaces; those with figural distances above the 92.5 percentile of the baseline distribution of figural distance are marked by ▴ and the rest, by ○. In B, we see that most of the marked episodes are caudal–lateral kicks of larger magnitudes. In D, we see that the marked episodes include both shorter rightward jumps and longer leftward jumps. These plots suggest that in selected regions of the behavioral space, there were load-related trajectory changes due either to incomplete compensation to the inertial load or to compensatory EMG changes that caused trajectory deviations. In B, the 2 boxes indicate the behavioral-space regions from which the kick examples shown in Figs. 4 and 5 are taken.

FIG. 4.

Examples of EMG and kinematic data of lateral kicks collected before, during, and after inertial loading. Shown here are 9 kicks with similar kick displacement vectors (their locations in the kick behavioral space shown in Fig. 3B) taken from preload (A), loaded (B–E), and postload (F–I) trials of the postload1, load2, and postload2 epochs of frog R, respectively. The temporal locations of these nine kicks within those epochs are indicated in the bar shown in the figure's top right corner. Each kick is divided into 2 phases, a and b, roughly corresponding to the extension and flexion phases, respectively. Phase boundary is indicated in each episode by a vertical dotted line. EMG data shown here (1,000 Hz) were high-pass filtered (finite impulse response filter, 50th order, cutoff of 50 Hz), to remove motion artifacts, and then rectified. Kinematic data shown here (59.64 Hz) were extracted from deinterlaced videoimages. Joint extension is indicated by an increasing angle (hip angle, dotted line; knee angle, solid line) and joint flexion is indicated by a decreasing angle. The EMG and kinematic data were synchronized by a digital counter. Here, the depicted EMG data are time shifted to the right by 62 ms to account for electromechanical delay. Figures in the bottom row display the ankle trajectory of each kick. In each of those bottom row figures, the origin denotes the initial ankle position; trajectory of phase a is shown in solid line and trajectory of phase b, in dotted line. Both the femur and tibiofibula are depicted as straight lines with lengths of 1. For more descriptions of these EMGs and kinematic data, see results (Trajectories before and after inertial loading and EMG changes after inertial loading).

FIG. 5.

Examples of EMG and kinematic data of caudal–lateral kicks collected before, during, and after inertial loading. Shown here are six kicks with similar kick displacement vectors (their locations in the kick behavioral space shown in Fig. 3B) taken from preload (A), loaded (B–D), and postload (E–F) trials of the postload2, load3, and postload3 epochs of frog R, respectively. The temporal locations of these 6 kicks within those epochs are indicated in the bar shown in the figure's top right corner. Specifications of the EMGs and kinematic data shown here are the same as those described in the legend to Fig. 4. Notice that the ankle trajectories of the unloaded kicks (A, E, F) rotate in the anticlockwise direction, whereas those of the loaded kicks (B–D) rotate in the clockwise direction. For more descriptions of these EMGs and kinematic data, see results (Trajectories before and after inertial loading and EMG changes after inertial loading).

FIG. 6.

Muscle synergies underlying the EMGs of kicking (A), jumping (B), stepping (C), and swimming (D) extracted by our NMF-based procedure. Each muscle synergy shown here, with vector magnitude normalized to unity, is represented as a bar graph across the 13 recorded muscles. The colors of the bars indicate the moment arm signs of the muscles (blue: hip extensor and/or knee flexor; red: knee extensor and hip flexor; black: major flexors in the thigh; green: ankle flexor and knee extensor; magenta: ankle extensor and knee flexor). Bars with component magnitudes >0.1 are labeled by their corresponding abbreviated muscle names (see methods, Implantation of EMG electrodes, for full names of the muscles). To the right of each synergy vector is a row of boxes for depicting a selection of loading conditions in which the synergy was found to be activated by the algorithm. Each box corresponds to a particular loading condition (as labeled in the topmost row in each figure panel). Those boxes corresponding to the conditions in which the synergy was found to be activated are colored gray and those boxes corresponding to the conditions in which the synergy was found to be inactivated are colored black. In A (kick) and C (step), similar pairs of synergies with complementary condition selections (i.e., their selections, when combined, span all conditions) are grouped together by a bracket. As can be seen, most synergies shown in the figure either have selections spanning all conditions or belong to a complementary group. These results suggest strongly that many muscle synergies were invariant across different unloaded and loaded conditions.

FIG. 7.

Reconstructing EMGs of kicking (A), jumping (B), and swimming (C) episodes using the extracted muscle synergies and their corresponding activation coefficients. In A, the kick episodes shown are the same as some of the kick examples shown previously in Figs. 4 and 5 (their corresponding figure labels shown in parentheses). The original EMG data (filtered, rectified, and integrated) are shown in thick black lines. Superimposed onto the EMG data are the reconstructed EMGs, decomposed into different colors denoting the contributions of the different synergies to the reconstruction, and matching the colors of the coefficient time traces shown below the reconstruction. The labels of the synergies (K1–K8, J1–J6, and Sw1–Sw7) also match those shown in Fig. 6. As described in results (Examples of EMG reconstructions), the plots here suggest that EMG changes associated with loading can be explained as altered activations of invariant muscle synergies.

FIG. 8.

Modulation of synergy activations after loading in selected regions of the behavioral space. A: before performing statistical tests for comparing of the unloaded and loaded coefficients, the onset and offset times of coefficient bursts were automatically marked by a burst detection procedure (see methods, Detecting activation bursts of synergy coefficients). Shown here is an example of a coefficient trace containing 5 bursts detected by our procedure. The onset and offset times found here agree very well with those marked manually using a Matlab graphical user interface. For each phase, statistical tests comparing unloaded and loaded burst attributes were performed on the maximum coefficient amplitude (in this example, p₂ for phase a and p₄ for phase b), the total burst duration (in this example, d₁ + d₂ for phase a and d₃ + d₄ + d₅ for phase b), and the burst onset times (phase b only; in this example, t₃, t₄, and t₅). B–F: in our comparisons of coefficient amplitude and burst duration, we compared the coefficients of loaded and unloaded episodes located at similar positions in the behavioral space. For each phase and for each baseline and loaded episode, the average difference between the amplitude or duration of that episode and those from its 5 closest baseline neighbors were calculated. The difference values of the loaded episodes were then compared against those from the baseline episodes and plotted in the behavioral space to assess whether modulation of the burst attribute was a function of the behavioral variety performed. In B to F, each data point represents a loaded episode; the shape of the marker indicates whether the difference value of an attribute for that episode lies below (▿), within (○), or above (▵) the interval containing 95% of the baseline difference values. In these 5 panels for different burst attributes, many loaded episodes have difference values lying above baseline. The triangles are also colored according to a blue-to-red scale to show their actual difference values. Our plots here show that a burst attribute (burst amplitude or burst duration) was modulated after loading only if the behavioral variety performed fell into selected regions of the behavioral space. See results (Modulation of synergy activations after loading in selected regions of behavioral space) for more details.

FIG. 9.

Inertial loading increases number of bursts and delays burst onset times. A and C: bar graphs showing the percentages of unloaded (gray) and loaded (black) episodes having particular numbers of coefficient bursts in phase b across the range of burst numbers observed. B and D: burst onset times of phase b bursts from unloaded (□) and loaded (○) episodes (mean ± SD). All burst onset times were specified as the time elapsed between the phase boundary and burst onset. Statistical significance of the difference between the unloaded and loaded onset times was assessed, for each order of burst, using the Kruskal–Wallis test (*P < 0.01).

FIG. A1.

An illustration of the steps involved in the NMF-based procedure we propose for discovering shared and specific features across an arbitrary number of data set conditions. In this example, there are data from 4 conditions, represented by a row of 4 square boxes. In each row, those boxes corresponding to the selected conditions (i.e., the conditions in which the synergy is activated) are colored gray and, if otherwise, colored black. Each block of boxes in the diagram shown under step 2 and step 4 represents the selections used in an extraction event; each row of the block denotes the condition selection of a synergy; and each column denotes the synergies selected for a condition. For descriptions and explanations of this figure, see the appendix (A procedure for discovering shared and specific synergies).

FIG. A2.

Validating the synergy extraction procedure by simulations. To test the performance of the procedure we developed, simulated EMG data (12 muscles, M1–M12) were generated by combining random synergies drawn from an exponential distribution (A). Data for 4 conditions were generated and each of the 7 random synergies shown here was assigned to a random condition selection (gray = on; black = off), generated in such a way that each condition comprises activations of 4 synergies. Our procedure was then applied to the simulated data. It returned both the synergies extracted and, for each, a condition selection specifying in which conditions the synergy is activated. The extraction results are shown here in B. Not only are the extracted synergies very similar to the original (scalar product = 0.9603–0.9968), but the extracted condition selections also perfectly match the original selections.

FIG. A3.

Quantification of the performance of our procedure in the simulated data sets. A: we computed the scalar product between best-matching pairs of original and extracted synergies as a measure of the procedure's performance. If a simulated synergy had more than one extracted synergy matched to it, its scalar product was averaged over the multiple matched pairs; if a simulated synergy was unmatched, a scalar product of 0 was assigned to it. We show here the scalar product distributions for both the unshuffled (thick solid line) and shuffled (dotted line) data sets. Whereas the scalar product values from unshuffled data all lie close to unity, those from shuffled data range from about 0.6 to 0.9. Clearly, our procedure performed very well in identifying synergy structures. We also extracted synergies from the pooled data set without applying our proposed procedure, initiating NMF with the total number of synergies in the simulated synergy set. Scalar product values from such extractions (thin solid line) overlap substantially with the baseline distribution, suggesting that our procedure enables better NMF performance. B: the procedure's ability to extract condition selections was assessed by counting the total number of selected and unselected conditions correctly identified for each simulated synergy. If multiple extracted synergies were matched to a single simulated synergy, the selections from the matched extracted synergies were combined; for the unmatched simulated synergies, the correctly identified number was set to 0. We see that most of the extracted selections from the unshuffled data sets (gray bars) contain 4 correctly identified selections (i.e., most extracted selections perfectly match their originals), whereas extracted selections from the shuffled data sets (black bars) spread over the entire range of numbers. Thus our procedure also performed very well in identifying condition selections of the synergies.