Abstract
Neuromuscular control of voluntary movement may be simplified using muscle synergies similar to those found using non-negative matrix factorization. We recently identified synergies in electromyography (EMG) recordings associated with both voluntary movement and movement evoked by high-frequency long-duration intracortical microstimulation applied to the forelimb representation of the primary motor cortex (M1). The goal of this study was to use stimulus-triggered averaging (StTA) of EMG activity to investigate the synergy profiles and weighting coefficients associated with poststimulus facilitation, as synergies may be hard-wired into elemental cortical output modules and revealed by StTA. We applied StTA at low (LOW, ∼15 μA) and high intensities (HIGH, ∼110 μA) to 247 cortical locations of the M1 forelimb region in two male rhesus macaques while recording the EMG of 24 forelimb muscles. Our results show that 10–11 synergies accounted for 90% of the variation in poststimulus EMG facilitation peaks from the LOW-intensity StTA dataset while only 4–5 synergies were needed for the HIGH-intensity dataset. Synergies were similar across monkeys and current intensities. Most synergy profiles strongly activated only one or two muscles; all joints were represented and most, but not all, joint directions of motion were represented. Cortical maps of the synergy weighting coefficients suggest only a weak organization. StTA of M1 resulted in highly diverse muscle activations, suggestive of the limiting condition of requiring a synergy for each muscle to account for the patterns observed.
SIGNIFICANCE STATEMENT Coordination of muscle activity and the neural origin of potential muscle synergies remains a fundamental question of neuroscience. We previously demonstrated that high-frequency long-duration intracortical microstimulation-evoked synergies were unrelated to voluntary movement synergies and were not clearly organized in the cortex. Here we present stimulus-triggered averaging facilitation-related muscle synergies, suggesting that when fundamental cortical output modules are activated, synergies approach the limit of single-muscle control. Thus, we conclude that if the CNS controls movement via linear synergies, those synergies are unlikely to be called from M1. This information is critical for understanding neural control of movement and the development of brain–machine interfaces.
Introduction
The manner by which the CNS coordinates the action of many muscles to control precise voluntary movements is a matter of debate. One simplifying mechanism potentially used by the CNS is the use of fundamental motor patterns or muscle synergies (Bernstein, 1967; Lee, 1984; Flash and Hogan, 1985; Hogan, 1985, 1988). A muscle synergy is a template of relative muscle activations used to create movement. Multiple muscle synergies could be activated simultaneously at varying levels, and each muscle could be represented in multiple synergy patterns. Synergies could be summed together to create complex movements (Ajiboye and Weir, 2009; Burkholder and van Antwerp, 2013; Steele et al., 2013). Previous studies have mathematically decomposed electromyography (EMG) data via non-negative matrix factorization (NMF) from a variety of voluntary movement tasks in humans and animals and found that the EMG of many (8–19) muscles can be represented by a relatively small number of fundamental muscle synergies (typically 2–5, but as many as 8 synergies; Jacobs and Macpherson, 1996; Saltiel et al., 2001; Ting and Macpherson, 2005; Klein Breteler et al., 2007; Overduin et al., 2008; Ajiboye and Weir, 2009; Torres-Oviedo and Ting, 2010; Dominici et al., 2011).
While NMF of voluntary EMG identifies potential synergies of voluntary movement, we applied an alternative approach using electrical stimulation to activate cortical pathways, which could elucidate how the CNS may be wired for synergies (Kasser and Cheney, 1985; Tresch et al., 1999; Saltiel et al., 2001; Mushahwar et al., 2004; Roh et al., 2011). High-frequency long-duration intracortical microstimulation (HFLD-ICMS) of the primary motor cortex (M1) was previously used to evoke EMG activity for synergy analysis (Amundsen et al., 2012; Overduin et al., 2012; Amundsen Huffmaster et al., 2017). We used single-electrode HFLD-ICMS to systematically explore the forelimb region of two rhesus macaques while recording EMG activity from 23–24 forelimb muscles (Van Acker et al., 2013, 2014; Amundsen Huffmaster et al., 2017). NMF analysis of this stimulus-evoked EMG yielded 4–5 synergies accounting for 90% of the variance in the EMG. Analysis of voluntary reaching data in the same experiment yielded 2–3 synergies. HFLD-ICMS facilitates muscle activations that cause movement of the forelimb, but this activation of multiple cortical networks around the electrode site may have obscured underlying fundamental output modules of M1. Since synergy results depend on EMG pattern variation (Ethier et al., 2006; Tresch et al., 2006), a more precise activation of M1 output modules using stimulus-triggered averaging of EMG activity (StTA; single-pulse ICMS; Cheney and Fetz, 1985; Kasser and Cheney, 1985; Fetz et al., 1989) allows for an investigation of potential synergies associated with the most fundamental organization of neural output networks within M1 (Cheney and Fetz, 1985).
The purpose of this study was to identify potential synergies and associated weighting coefficients present in putative fundamental facilitatory output modules of M1. Synergies and weighting coefficients were extracted using NMF of the mean EMG poststimulus facilitation (PStF) evoked by StTA applied systematically throughout the entire forelimb representation of M1. Poststimulus suppression was not included in the analysis. StTAs were taken at the same cortical locations previously studied with HFLD-ICMS (Van Acker et al., 2014; Amundsen Huffmaster et al., 2017). StTA data were collected from two rhesus macaques performing a reach-and-prehension task. Given the expectation of a much larger range of EMG activation patterns with StTA compared with HFLD-ICMS or voluntary movement, we hypothesized that the number of NMF synergies needed to account for the StTA activation patterns would be greater and might approach single-muscle control, that is, one synergy per muscle. While single-muscle control was not quite necessary, the characteristics of the synergies and their weighting coefficient maps suggested a trend toward nearly singular muscle activation patterns.
Materials and Methods
Experimental design and statistical analysis.
Single microelectrode StTA was systematically applied throughout the right forelimb region of (left brain) M1 in two awake male rhesus macaques during a reach-and-prehension task. Trains of electrical stimuli were applied to each of 247 putative lamina-V sites in M1 (140 sites in Monkey A; 107 sites in Monkey X) at LOW (∼15 μA) and HIGH (∼110 μA) intensities. Corresponding EMG was recorded from 24 forelimb muscles (Table 1) and StTAs were calculated during data collection and again post hoc for each cortical location. NMF was applied to the mean PStF of the StTA data, with a single StTA trial defined as a single average of hundreds of stimulus pulses. Synergies were selected for each of the four datasets (Table 2) using a 90% variation-accounted-for (VAF) threshold. Synergy profiles were compared using Pearson's correlation coefficients. The weighting coefficients were mapped to the electrode location in M1; the centers of the maps were matched across monkeys and the resulting matrix of weighting coefficients was compared using two-dimensional correlations.
Muscles recorded in EMG, listed by their primary joint of action
Dataset characteristics
Subject preparation.
During training and data collection, each monkey was seated in a custom chair designed to comfortably restrain the left forelimb at his side. A waist plate was positioned to keep the hindlimbs out of the testing space. A clear acrylic shield was placed in front of each monkey's face with a small hole to allow him to feed himself (Park et al., 2004). Each monkey was trained to perform a right-forelimb reach-and-prehension task, which consisted of the following: (1) depressing a home plate lever at waist level for 1 s to activate a food reward, (2) retrieving a banana pellet by extending the hand to a food well at shoulder level, (3) bringing the pellet to the mouth, and (4) returning to the home plate (Park et al., 2004; Van Acker et al., 2013).
Between training and data collection, each macaque underwent surgical implantation of a titanium recording chamber over the left hemisphere. The chamber was positioned stereotaxically, centered over the forelimb region of M1, and attached with orthopedic titanium bone screws. An X-Y positioner was secured to the chamber and used to lower single glass-insulated platinum-iridium electrodes (FHC; impedance, 0.7–1.5 MΩ) through the dura and into lamina V of the cortex for recording neural activity and performing stimulation (McKiernan et al., 1998, 2000; Park et al., 2001).
EMG electrode wires were implanted subcutaneously in 24 forelimb muscles, extending from a cranial connector (Amphenol). Two wires, spaced ∼5 mm apart, were inserted to a depth of ∼2–3 cm into the belly of each muscle. Placement was confirmed via muscle twitches observed from short trains of stimulation (Natus Medical, Grass SD9 stimulator; McKiernan et al., 1998, 2000; Park et al., 2001; Hudson et al., 2010). Cross talk was tested by compiling EMG-triggered averages of EMG activity from the trigger muscle and all nontrigger muscles (Buys et al., 1986; Van Acker et al., 2013). All muscles had negligible cross talk except muscle 24 (first dorsal interosseous muscle) in Monkey X, which was excluded from the database. All surgeries were performed as directed by the standards of the Guide for the Care and Use of Laboratory Animals, published by the U.S. Department of Health and Human Services and the National Institutes of Health, and the protocol was approved by the Institutional Animal Care and Use Committee.
Data acquisition.
While the macaque was engaged in the reaching task, EMG data were recorded from 24 forelimb muscles and trains of electrical stimuli were applied to lamina-V sites in M1. Electrode tracts were spaced at 1 mm intervals along the surface of the precentral gyrus and at 0.5 mm intervals in tracks extending vertically down the precentral gyrus (Cheney and Fetz, 1985; Van Acker et al., 2013). Stimulation pulses were biphasic square waves consisting of an initial negative pulse of 0.2 ms followed by a positive 0.2 ms pulse. EMG data were filtered from 30 Hz to 1 kHz and digitized at 4 kHz (Cambridge Electronic Design, Power 1401). Placement of the electrode in lamina V was based on the presence of large neuronal spikes as well as confirmation of visible PStF in StTAs at 15 μA and 10 Hz (Park et al., 2001). If PStF was not observed, the electrode was raised or lowered until effects were obtained. Data from locations lacking PStF were discarded. Some sites presented difficulties in identifying layer V, so the cortical maps (see Figs. 5, 6) show a few empty (gray) boxes and fewer tracks in the sulcus for Monkey X. Tracks in which layer V remained elusive were reattempted on a later day.
At each cortical site, StTA data were collected at 10 Hz, first using low (LOW; typically 15 μA) intensities to represent the most focal and minimal activation of the neural network. LOW StTA was collected first to help confirm a layer-V site for further data collection. In a previous study we showed that the effects of stimulation are not dependent on stimulation order (Griffin et al., 2009). These responses were compared with a high-intensity StTA collected next (HIGH; typically 110 μA). The StTA data for this analysis and the HFLD-ICMS data published previously (Amundsen Huffmaster et al., 2017) were collected as part of the same protocol. The HIGH StTA intensity was similar to the intensity used for HFLD-ICMS and was collected immediately following the LOW StTAs and before the HFLD-ICMS data collection. Because the effects from low-intensity StTA are generally weaker, more epochs (stimuli) are usually required to clearly resolve the baseline from the poststimulus effects. On average, ∼2000 stimulus pulses were used for LOW StTA and 1000 for HIGH StTA. The presence of PStF was assessed in near real time during data collection on a separate computer and recalculated post hoc for final analysis. Thus, StTA data were collected until poststimulus effects were clearly identified or until it was determined that we were not at a layer-V site, at which point a different depth was tested.
EMG processing and StTA analysis.
The raw EMG data were loaded into Matlab (Mathworks), converted from integers to voltage values, and full-wave rectified about the mean. Stimulus artifacts were removed by replacing artifact data with baseline data for 2 ms poststimulus.
StTAs were calculated for each muscle at each cortical location from a long train of stimulus pulses (average, ∼1500 stimuli). The data train was sorted into epochs of 60 ms of EMG data surrounding each stimulus pulse (20 ms before the trigger to 40 ms after). Each epoch therefore had 240 data points, corresponding to a bin width of 250 μs, at the collection frequency of 4000 Hz (Cheney and Fetz, 1985; Park et al., 2001). The epochs for a trial were averaged over each individual 60 ms epoch. An epoch was included in the StTA when the muscle's average EMG in that epoch was ≥4% of the full-scale signal (McKiernan et al., 1998; Park et al., 2004). Very few epochs were discarded based on this criterion.
Onset and offset of PStF effects were identified when the StTA crossed the pretrigger baseline mean by ±3 SDs for ≥6 data points (1.5 ms; Park et al., 2004). Pure poststimulus suppression (PStS) effects were identified in a similar manner. Only the facilitation component of biphasic effects (facilitation followed by a trough) was counted as a PStF effect because of ambiguity concerning the latency and origin of the trough (Park et al., 2004). The magnitude of PStF was calculated for each muscle as the mean percentage facilitation between onset and offset, using the following equation (Cheney and Fetz, 1985; Kasser and Cheney, 1985):
This is a widely used measure of the strength of effects on different motoneuron pools (Cheney and Fetz, 1985; Kasser and Cheney, 1985; Buys et al., 1986; Baker and Lemon, 1998; Schieber and Rivlis, 2005; Davidson et al., 2007). The magnitude of PStF is independent of baseline level and is highly stable across a wide range of behavioral and experimental conditions (Griffin et al., 2009). One explanation for this is that cortical cells tend to make synaptic contact with all the motoneurons in a motoneuron pool (Palmer and Fetz, 1985). As more motor units are added to a movement in association with greater voluntary effort, the PStF signal grows approximately in proportion to the baseline level of EMG activity (Fourment et al., 1995). Taking all this into account, StTA is highly effective in identifying the targets and strength of corticospinal connections to motoneuron pools. StTA synergies reveal how cortical sites facilitate combinations of muscles as potential synergies.
Another issue concerns the possibility of multiple deflections or effects in the same average for one muscle. For instance, there could also be a delayed secondary facilitation peak after the initial facilitation peak. In these cases, which are rare, we used the maximum facilitation peak, which was almost always the first peak, and we ignored the secondary facilitation peak. Also, biphasic effects, consisting of initial facilitation followed by a trough, are common. Troughs following a facilitation peak were not counted as suppression because of ambiguity about the origin of the trough. Although most are probably real suppression effects coming in at a longer latency than facilitation, such troughs could also be related to the refractory period of the motor unit.
One-way ANOVAs were used to investigate any differences between the percentages of poststimulus facilitation effects in each muscle grouping (groups at the shoulder, elbow, and wrist/digits/hand). Bonferonni corrections were used for multiple comparisons, with a significance value set to p = 0.05.
NMF.
NMF was used to determine the synergistic patterns of muscles and cortical organization of these synergies as activated by localized stimulation of M1. Because NMF depends on input without negative values, PStS effects were identified, but set to zero, and NMF was performed on each dataset of mean percentage facilitation from the LOW-intensity and HIGH-intensity datasets for each monkey. Pure PStS effects represented 30.7% of the total effects in our LOW dataset, and 15.9% of the total in the HIGH dataset. While inhibitory output from the cortex to motoneurons is very important, we chose to use NMF to preserve the ability to examine our results within the context of our previous study and the studies of others. NMF uses an iterative update method to factor a matrix of mean percentage facilitation values into two matrices representing synergies and weighting coefficients (Paatero and Tapper, 1994; Lee and Seung, 1999, 2001; Torres-Oviedo and Ting, 2010). Synergies are vectors containing non-negative values for each muscle denoting the likelihood that groups of muscles will be activated simultaneously at consistent ratios relative to one another. These vectors are nonorthogonal and normalized to have one muscle activated to a level of one. The weighting coefficients are vectors containing a value specifying the level at which a synergy was activated for each trial in the original data. Together, the EMG pattern is represented as follows: EMG = V = W * H (Lee and Seung, 2001) where V is a (m × t) matrix of original (or later reconstructed) mean percentage facilitation values, with t trials of StTA and m muscles, W is a (m × k) matrix of k column synergy vectors, and H is a (k × t) matrix of weighting coefficients with one value showing how each trial used each synergy (Lee and Seung, 1999; Donoho and Stodden, 2003; Berry et al., 2007; Torres-Oviedo and Ting, 2007). NMF was used to factor the original dataset into sets of synergies and weighting coefficients using the multiplicative update rule (Paatero and Tapper, 1994; Lee and Seung, 1999, 2001; Tresch et al., 2006; Berry et al., 2007; Ajiboye and Weir, 2009; Badeau et al., 2010; Burkholder and van Antwerp, 2013).
Rank selection: VAF.
NMF was repeated for every potential number of synergies, or rank, from 2 to 22. The minimum number of synergies needed to describe ≥90% of the total VAF from the original data was used for further analysis. VAF is the uncentered Pearson correlation coefficient: VAF = 1 − (SSE/TSS), where SSE is the sum of the squared errors and TSS is the total sum of the squares, taken with respect to zero (d'Avella et al., 2006; Torres-Oviedo et al., 2006; Torres-Oviedo and Ting, 2010; Amundsen Huffmaster et al., 2017). Synergies are presented from the algorithm in an unweighted and randomized order and the VAF of the whole synergy set is given. VAF is not available for each individual synergy.
Synergy profiles and weighting coefficients.
Synergy vectors were normalized to the maximum muscle activation in each synergy and plotted as bar graphs and radial plots, with the muscles ordered according to their actions. When comparing sets of synergies, synergy vectors were matched and sorted in the order of best match by using the minimum dot product between any two synergies. Once the synergies were matched, they were compared using correlation coefficients and the corresponding p values, with significance set at p < 0.05 (Overduin et al., 2012; Amundsen Huffmaster et al., 2017).
The weighting coefficient set for each synergy was mapped to the corresponding M1 electrode location by plotting a color-coded square at each site to represent how strongly a particular synergy was weighted to explain the StTA results obtained at that site. We attempted to compare the synergy maps between the two monkeys. To do this we flattened the maps to two dimensions for simplicity, aligned the maps from each monkey at the central point of the sulcus, and matched each set of two points counting radially outwards from that central point (see Figs. 5, 6, white asterisks). Due to the different shapes in the maps, we were able to create an 8 × 8 grid of 64 matched site locations (45% of Monkey A's sites and 59% of Monkey X's sites). We calculated the two-dimensional correlation coefficient (corr2 in Matlab) for our matched grid to compare the two monkeys' maps for each synergy weighting coefficient set, including sites within the white boundaries in Figures 5 and 6. To estimate a significance value, we used a resampling technique to test the null hypothesis that the maps were more similar than would occur randomly. To do this, we used permutation testing where we randomly assigned the values of each map to a new site and calculated the two-dimensional correlation coefficient across the two maps. We repeated this process 10,000 times to account for the Monte-Carlo variance, and we estimated a p value as the proportion of comparisons that were greater than the calculated value using the real maps. Significance was set to p < 0.05 (Josse et al., 2008).
Results
Dataset and synergy characteristics
The general characteristics of the four datasets are shown in Table 2. Within any trial (a single average of hundreds of stimulus pulse responses), PStF occurred on average in 3.8 (Monkey X) and 5.9 (Monkey A) muscles for LOW trials and 11.3 (Monkey X) and 14.6 (Monkey A) for HIGH trials. An example of StTA applied to a single site is shown in Figure 1. Low-intensity StTA elicited fewer poststimulus effects than the high-intensity StTA. One-way ANOVAs on the percentage of PStF in each muscle group (shoulder, elbow, and wrist/digits/hand) showed a main effect of functional muscle group for the LOW-intensity StTA datasets in both monkeys (A: F = 3.5, p = 0.047; X: F = 4.7, p = 0.022), but no effect for the HIGH datasets (A: F = 0.49, p = 0.621; X: F = 1.7, p = 0.200). Follow-up comparisons using the Bonferroni correction for multiple comparisons showed no significant differences between muscle groups for A LOW, but X LOW showed that the shoulder muscles were significantly less likely than the wrist/digits/hand muscles to have a PStF (p = 0.022).
EMG responses at a single cortical site derived from StTA. StTAs were collected using a stimulus intensity of 110 μA and a frequency of 10 Hz, and were applied throughout the reaching task. The stimulus event onsets are denoted by the thin vertical line. StTAs were based on ∼ 1000 trigger events and were processed to remove stimulus artifacts. For muscle abbreviations, see Table 1.
Synergies extracted using NMF on the PStF of StTA showed that 10 or 11 synergies were required to account for ≥90% of the original LOW StTA data for Monkeys A and X, but only 4 or 5 synergies were needed for the HIGH StTA dataset (Figs. 2, 3; Table 2).
VAF with increasing number of synergies for the two StTA datasets for each monkey. The 90% threshold is shown, and the first synergy rank to cross above that threshold is highlighted by a red box.
A, B, Muscle synergies are shown for the two StTA datasets for (A) Monkey A and (B) Monkey X. The LOW StTA was compared with the HIGH StTA synergies with the correlation coefficients (R) written beneath the plots. * denotes p < 0.05. Synergies are arranged in order of the best to worst match. Synergies without a match are presented in random order, as they are given from the NMF algorithm. Each muscle is represented by a single bar in the order listed at the right (the first dorsal interosseous muscle was only present in Monkey A).
Synergies compared: low-intensity versus high-intensity StTA
All HIGH synergies except one for each monkey were significantly correlated (p < 0.001, Pearson's correlation) to a synergy in LOW-intensity datasets (see Fig. 3 for R values). The mean correlation coefficient between synergies was R = 0.71 ± 0.4 (mean p = 0.10) for Monkey A and R = 0.65 ± 0.45 (mean p = 0.25) for Monkey X. In general, the LOW synergies showed a sparse muscle representation, with only one or two muscles at middle to high activation within each synergy, while the HIGH synergies generally had substantially more muscles (1–6) activated at middle and high levels.
Synergies compared across monkeys
The synergy profiles were compared across monkeys and plotted with the muscles in each synergy arranged radially around a unit circle (Fig. 4). Of the LOW synergies (10 for A; 11 for X), nine synergies were significantly correlated (p < 0.01) between monkeys. When comparing the synergies from the high-intensity datasets across monkeys (4 for A; 5 for X), three synergies were significantly correlated (p < 0.01).
A, B, Radial plots of StTA synergies compared across monkeys for (A) LOW and (B) HIGH current intensity datasets. Synergies are from the top two rows of the bar plots in Figure 4. C, The legend shows how the muscles are organized around a unit circle with muscles arranged in 15° increments. The outermost circle has value of 1, the middle circle is 0.5, and the two dashed circles show a 0.25 and 0.75 level of muscle activation. Muscles with a high activation (>0.5) in a synergy are labeled on each plot. Pearson correlation coefficients (R) and corresponding p values show the relationships between the two monkeys. Synergies are arranged in order of the best to worst match across monkeys for each of the LOW or HIGH StTA sets.
Muscle and movement representation
Most of the LOW synergies (5–6 of 10–11) and even some of the HIGH synergies (1–2 of 4–5) consisted of a single muscle represented at a high level while several other synergies consisted of a few muscles acting as close synergists. In the LOW synergies from Figure 4, single-muscle representation included synergies 1 (flexor digitorum profondus), 2 (dorsal epitrochlearis), 3 (triceps lateral head), 5 (extensor carpi radialis in X), 10 (brachialis in A, flexor carpi radialis in X), and 11 (teres major in X). In the HIGH dataset, synergies 1 (flexor digitorum profondus) and 3 in A (dorsal epitrochlearis) had single-muscle representation. Synergies involving only a few muscles in the LOW group were synergy 4 (finger extensors), 5 (wrist extensors), 7 (finger and wrist extensors), and 9 (finger and wrist flexors), and in the HIGH dataset were 2 (finger extensors) and 3 (elbow extensors). In contrast, only a few synergies had representation of disparate muscles: in the LOW group these were only synergies 6 and 9 (in A only) and in the HIGH dataset were synergies 4 and 5 (in A only).
Most of the LOW synergies had strong representation of the hand/wrist/digits and/or the elbow. Despite the similar presence of PStFs in different muscle groups for Monkey A (Table 2), only Monkey X had a single synergy with strong shoulder activity (teres major, which adducts and medially rotates the arm). While most synergies in the low-intensity set showed single-muscle or nearly single-muscle representation, two of Monkey A's synergies had representation in both the elbow and a distal joint. The muscles represented a wide range of directions, except for the notably absent shoulder activity for most directions in both monkeys.
Of the HIGH synergies, synergy 1 primarily flexed the digits, synergy 2 extended the elbow, and synergy 3 extended the digits. Strong multijoint activations were only seen for Monkey A in synergy 5 with flexion of the elbow and digits. Synergy 4 in Monkey X activated both flexors and extensors in the hand and digits. Synergy 4 for Monkey A was primarily wrist extension. Similar to the synergies from low intensities, the synergies represented a wide range of directions for movement of the lower arm, but there were no shoulder muscles represented in the HIGH synergies (only one shoulder muscle was represented in the LOW synergies).
Weighting coefficient characteristics
The weighting coefficients for synergies from each of the low-intensity and high-intensity StTA datasets (Fig. 4, radial plots) were mapped to the stimulus location in the brain (Figs. 5, 6) with color-coded squares representing the amount each synergy was weighted to explain StTA results at that location. The similarities between the two monkeys' maps (based on a common central zone of each map) were computed as correlation coefficients between the matched regions. Significance values (p < 0.05) were estimated via permutation testing of the two-dimensional correlation coefficient. The actual correlation coefficients and bootstrapped p values comparing monkey maps for each synergy are shown in Figures 5 and 6.
Mapping coefficients for synergies from LOW-intensity StTA. Heat maps showing the average weighting coefficients for each synergy at each stimulation location. Synergies are ordered to match the radial plots (arranged in order of the best to worst matched synergies across monkeys). The maps were compared across monkeys using two-dimensional correlations. As the anatomies of the two monkeys differ, the maps were matched at the center (white asterisk on the maps) and then the sites were aligned in a matrix moving outwards from that star. Sites included in the two-dimensional correlations are inside the white boundary lines. R and estimated p values are shown to compare the similarity between monkeys for each set of maps.
Mapping coefficients for synergies from HIGH-intensity StTA. See the Figure 5 legend for more details.
The LOW set for each monkey contained some maps with only one or two sites with high weighting coefficients (Fig. 5, A1, X1), while other maps had a large range of weights represented across the entire mapped region (A10, X10). The LOW synergy profiles 1–9 were found to match closely between the two monkeys and, while the weighting coefficient Maps 1–4 and 7–10 also appear to match qualitatively across monkeys, only Maps 1, 2, and 8 were more significantly correlated than random (p < 0.05). Most of the maps that matched across monkeys had localized points with very high weighting activation at locations very near each other (i.e., LOW Maps A1 and X1).
To check for uniform distributions in our permutation testing, histograms (data not shown) were created from randomly shuffling of the points of the map of Monkey A and then recalculating the R value to the map of Monkey X. This was repeated 10,000 times. The p value was then estimated from the proportion of the shuffled R values that were greater than the real R value. If the real R value was >5% of the shuffled-data R values, then it was declared significant. However, when inspecting histograms of the R values generated by the 10,000 permutations (standard in the field, data not shown), 8 of our 10 histograms for the LOW maps appeared intensely right-skewed, which could have artificially inflated our significance. We calculated the Fisher–Pearson third moment coefficient of skewness (Sk), adjusted for sample size, where the normal distribution has an Sk = 0. We also calculated the kurtosis (K) as the fourth central moment of X, divided by the fourth power of its SD, adjusted for bias, so that a normal distribution would have a K = 3 (Doane and Seward, 2011). Our histograms for the LOW dataset had Sk values ranging from 0.4 to 2.5, with K values ranging from 2.9 to 10.9, making interpretation of our permutation testing more challenging. Of the three significantly correlated LOW maps (Maps 1, 2, and 8), Maps 1 and 8 were right-skewed (Sk = 2.5 and 1.9; K = 10.9 and 6.8), while Map 2 was more normally distributed (Sk = 0.4; K = 2.9). The lack of normal distribution in general is probably due to the sparsity of the data matrix on the maps, as very few sites had high activations.
The weighting coefficient maps for the HIGH StTA synergies followed similar trends as the LOW maps, but in general showed larger clusters of high activation of each synergy (Fig. 6). Most of the maps had localized clusters, but HIGH X4 had a rather wide dispersion of high activations, including the sulcus and the gyrus. The fact that this synergy included relatively strong activations of several distal muscles (Fig. 4) undoubtedly contributed to this broad representation. Maps 1–3 were significantly more similar across monkeys than random (p < 0.05), but maps for A4 and X4 were quite different. The histograms for the HIGH Maps 1 and 2 appeared fairly normally distributed (Sk = 0.5 and 0.4; K = 2.9 and 2.3), and Maps 3 and 4 showed histograms appearing to be very close to normally distributed (Sk = 0.2 and 0.2; K = 2.8 and 2.6), giving more confidence in the permutation testing for the HIGH maps than the LOW.
The spatial organization of the weighting coefficient maps was qualitatively similar to that of maps from previous StTA mapping studies (using both different and the same monkey subjects) when including muscles similar to those represented in each synergy (Park et al., 2001; Van Acker et al., 2013). In both LOW and HIGH maps, the synergies activated in the sulcus consisted of more distal muscles, while those activated on the surface of the gyrus consisted of proximal or proximal plus distal muscles. Most of the LOW maps show one or two sites of stimulation that were very high with little gradation to the surrounding sites that had low to no activation, while HIGH maps tended toward better gradation between sites of high activation and those of low activation. The third type of map exhibited in Figures 5 and 6 had a wide range of activation across the forelimb region of M1 with no discernible pattern, as in LOW maps A10 and X10. These maps had high and mid-level activations throughout both the surface of the precentral gyrus and the bank of the central sulcus. LOW synergy 10 for Monkey A activated the brachialis, an elbow flexor, and LOW synergy 10 for Monkey X activated the flexor carpi radialis, a flexor of the hand.
Discussion
We investigated muscle synergies represented in output modules throughout the full extent of the M1 forelimb representation in two rhesus macaques as calculated via NMF applied to PStF from LOW-intensity (∼15 μA) and HIGH-intensity (∼110 μA) StTA. Muscle synergies could offer simplified neuromuscular control (Bernstein, 1967; Overduin et al., 2012), but the neural origin of synergies remains elusive (Giszter et al., 1993; Tresch and Bizzi, 1999; Saltiel et al., 2001; Mushahwar et al., 2004).
NMF from PStF required 10–11 synergies to account for 90% of the variance in the LOW StTA and 4–5 for the HIGH StTA. Synergies from LOW StTAs typically activated one or a few close synergistic muscles. While HIGH StTA synergies included more muscles, some consisted of a single muscle or a few synergists. Some synergies were consistent across monkeys, while others were subject-specific. All joints were represented, but not all the directions of motion could be achieved with the synergies. The lack of shoulder representation in the LOW dataset might be explained by the relative lack of PStF in shoulder muscles. However, this was not the case in the HIGH dataset. It is possible that the shoulder representation in the HIGH dataset synergies were part of the 10% variance left unaccounted for in each synergy set. Cortical maps of the weighting coefficients showed localized activation patterns for the synergies. StTA may be activating a synergistic network within M1 or further downstream, but the synergy profiles and weighting maps cannot confirm that these synergies represent a simplified control system.
Sparse muscle synergies with StTA
There are important differences between StTA-evoked EMG, HFLD-ICMS-evoked EMG (from previous work), and voluntary movement-evoked EMG (Amundsen Huffmaster et al., 2017). StTA causes no overt movement of the arm and generally no observable EMG response without averaging. Each single pulse of stimulation will excite corticospinal neurons at the site of stimulation and probably some interconnected neurons at a distance (Slovin et al., 2003; Tehovnik et al., 2006). Low-intensity StTA minimizes activation of the neural network, allowing for higher resolution of cortical output organization (Omrani et al., 2017). On the other hand, the high frequency and higher intensity of HFLD-ICMS produces much broader activation of the neural network overriding (hijacking) natural cortical activity and causing overt stimulus-specific movement of the forelimb in space (Griffin et al., 2014). Voluntary movement-related EMG is, by definition, driven by natural cortical network activity.
Direct comparison of StTA data to HFLD-ICMS data is complicated by differences in EMG normalization procedures. PStF was measured as a percentage over baseline, which removes the underlying level of muscle activity related to task performance, while the EMG was directly measured and normalized by session in the HFLD-ICMS and voluntary movement protocols. Given this complication, we have avoided direct comparison of StTA results with those from our previous work (Amundsen Huffmaster et al., 2017).
However, observationally, as the level of neuronal activation increased from LOW StTA to HIGH StTA to HFLD-ICMS (Amundsen Huffmaster et al., 2017), the sparsity of the muscle representation in each synergy set decreased and the region of synergy activations in M1 increased. The HFLD-ICMS and voluntary synergies often had activations of muscles across two or three joints, with every joint direction represented in ≥1 synergy. HFLD-ICMS-weighting maps had larger regions of cortical representation compared with StTA-weighting maps, but did not seem as functionally organized as other studies observed (Overduin et al., 2012, 2014, 2015). As cortical output resolution was augmented using more focal activation methods (i.e., LOW StTA), our synergy analysis approached the limit of requiring a synergy for each muscle to account for the patterns observed. In the case of such sparse synergies, any NMF-derived synergies would lose the advantage of simplifying movement control, even if these synergies did exist within a neural representation.
Synergy origin
While simplified control schemes are useful for engineering design, complex control neural networks may have developed habitually. The CNS could use a locally optimal control strategy to create habits for regularly used movements. Habitual synergies would evolve gradually (de Rugy et al., 2012) as skilled movements developed (von Hofsten and Rönnqvist, 1988; Myer et al., 2005; Hug et al., 2010; Dominici et al., 2011), in response to the individual's environment and experiences in the way that motor map topography depends on the subject's particular movement performance (Nudo et al., 1992, 1996; Graziano et al., 2004). They would not be highly organized or accessible through stimulation (Andersen et al., 1975). Control would be subject-specific and only similar due to shared experiences or physiology.
Habitual versus optimal control may be considered in how the cortex is organized for control of functional groups of muscles in two complementary and opposite ways: multiple corticomotoneuronal cells project to a single motoneuron (convergence) and a single corticomotoneuronal cell projects to multiple motoneurons (divergence; Andersen et al., 1975; Shinoda et al., 1981, 1986; Cheney and Fetz, 1985; Li and Martin, 2002). Spike-triggered averaging, similar to StTA, identifies the target muscles of single corticomotoneuronal cells, viewed as the smallest unit of cortical output (Fetz and Cheney, 1980; Lemon et al., 1986; Schieber, 1995). Simple synergies, such as cofacilitation of agonist muscles and reciprocal suppression of antagonist muscles, are hard-wired into corticospinal output networks (Kasser and Cheney, 1985). However, the hand region of M1 (Penfield and Boldrey, 1937) has overlapping neuronal populations for each digit and movement direction (Schieber and Hibbard, 1993). Complicated control should be expected from biological creatures with inherent redundancy and inefficiencies.
It seems unlikely that NMF-derived synergies are hard-wired into the neural networks of M1, as our cortical maps did not show localizations of the synergies. M1 synergies may be partial representations of motor commands activated elsewhere, such as the premotor area activating multiple M1 areas (Woolsey et al., 1952), the cerebellum initiating compound movements (Overduin et al., 2012; Thach, 2014), or motor outputs combining in the spinal cord (Giszter et al., 1993; Hart and Giszter, 2010; Levine et al., 2014; Caggiano et al., 2016). Synergies may also be encoded in the brainstem (Roh et al., 2011) or spinal cord (d'Avella et al., 2003; Overduin et al., 2009; Desmurget et al., 2013; Van Acker et al., 2013; Giszter, 2015; Graziano, 2016; Saltiel et al., 2016; Takei et al., 2017).
Study limitations
This study must be considered within the limitations of NMF and the study protocol (Amundsen Huffmaster et al., 2017). NMF assumes that synergies exist in the data in linear combinations with only positive values. This may or may not be reasonable (Paatero and Tapper, 1994; Lee and Seung, 1999; Ethier et al., 2006; Ajiboye and Weir, 2009; Takei et al., 2017), and is why we set any PStS to zero. NMF is dependent on the number and choice of muscles tested (Steele et al., 2013), the constraints the limb encounters (Burkholder and van Antwerp, 2013), and the variety of movement tasks and EMG patterns measured (Ethier et al., 2006; Tresch et al., 2006). We measured more muscles than is typical and used StTA to activate small clusters of corticomotoneuronal cells (Cheney et al., 1985) to increase output diversity (Torres-Oviedo and Ting, 2007; Dominici et al., 2011; Overduin et al., 2012; Steele et al., 2013, 2015). Cortical mapping was more comprehensive than done previously (Overduin et al., 2012, 2014, 2015). While maps were sometimes correlated, only central sites were included in the correlation analysis and the sparsity of data on the LOW maps make drawing conclusions challenging. Instead of applying StTA at each site's motor threshold, we used a consistent intensity to reduce variables and time. Alternative procedures, including algorithms that allow input of negative suppression values, are left to future studies.
Conclusion
Our goal was to investigate the presence of muscle synergies in M1 using StTA facilitation of EMG activity. The synergy profiles and weighting coefficient maps are difficult to interpret as a complete and simplified control scheme for forelimb movement. StTA produced very diverse muscle activations with a clear trend of increased synergy numbers with increased output diversity. LOW StTA required more synergies comprised of one or a few muscles, trending toward the limit of a synergy for each muscle and a loss of the simplified control scheme. This adds to a body of evidence raising doubts about the physiological significance of synergies derived from data-reduction approaches (Mushahwar et al., 2004; Kirsch et al., 2014; Mollazadeh et al., 2014; Liu et al., 2015). M1 motor control may not be optimal; it could be habitual, synergistic from elsewhere in the CNS, or a combination of the two. In M1, the CNS appears to select small assemblies of neighboring neurons with similar patterns of connections to motoneuron pools based on a match between the connections of these neurons and knowledge of the movement goal (Cheney and Fetz, 1985). These neurons represent the fundamental units of cortical output specifying hard-wired functional muscle synergies, which the CNS must select to produce the desired movements (Omrani et al., 2017).
Footnotes
This work was supported by National Institutes of Health (NIH) Grants NS-051825 and NS-064054, NIH Center Grant HD-02528, the University of Kansas Medical Center Biomedical Research Training Program, and the University of Kansas Madison and Lila Self Graduate Fellowship. We thank Ian Edwards, University of Kansas Medical Center; Molly McVey, PhD, University of Kansas Mechanical Engineering Department; and Annaria Barnds, PhD, University of Kansas Bioengineering Department for their technical assistance, along with Dr. Colum D. MacKinnon and the rest of the Movement Disorders Laboratory team at the University of Minnesota for their support. We also acknowledge funding from S. Amundsen Huffmaster's current institution, the University of Minnesota, including the MnDRIVE Postdoctoral Fellowship, NIH Grants R01 NS085188 and R01 NS088679-02, and the UMN Udall Center: National Institute of Neurological Disorders and Stroke Grant 1P50NS098573-01.
The authors declare no competing financial interests.
- Correspondence should be addressed to Dr. Paul D. Cheney, Department of Molecular and Integrative Physiology, University of Kansas Medical Center, Kansas City, KS 66160-7336. pcheney{at}kumc.edu
