For a given *N*I or *N*G, the algorithms iteratively updated the structures *ν*I *n* or *ν*G *n* and coefficients *w*I *n* and *w*G *n* until the total reconstruction error, *R*2, increased by less than 0.001 over 10 iterations. The algorithms were repeated five times for each extraction; the set of synergies with the highest EMG variation explained was selected for further analysis. Dimensionalities (*N*I and *N*G) were chosen by applying a threshold of *R*2 = 95% (Overduin et al., 2012).

### **MATCHING**

Synergies *ν*I *n* were compared and matched to synergies *ν*G *n* using a greedy search procedure (Tresch et al., 1999; Overduin et al., 2012). For all *N*I × *N*G possible pairs of ICMS-evoked and grasping-related synergies, we first computed dot products (e.g., 6 × 8 = 48 dot products, for G2. The pair of ICMS-evoked vs. grasping-related synergies with the highest dot product was defined as the best-matching pair. The pair with the highest dot product among the remaining (*N*I − 1) × (*N*G − 1) pairs (e.g., 5 × 7 = 35 pairs, for G2) defined the second-best match. This process was repeated until all synergies in one set had been paired (e.g., over min(6,8) = 6 times for G2). We then used Monte Carlo simulation to assess the significance of each match. We repeated the greedy search algorithm 10,000 times for each monkey, after first randomly shuffling EMG channel identity each time. (For G2, for instance, this involved finding 10,000 × 48 = 480,000 dot products.) We then compared the highest (best-matching) dot product between *actual* ICMS-evoked and grasping-related synergies with the distribution of highest dot products from the 10,000 comparisons of *shuffled* synergies. If the former value exceeded the 95th percentile of the distribution of latter values, we took the match as significant at *p* < 0.05. We then repeated this comparison for the second-best actual synergy pair vs. the distribution of second-best shuffled pairs, etc.

averaged over ICMS trains.) Of G2's ICMS sites, this synergy was most strongly evoked at site XI, in MI.

### **ANALYSIS**

In reconstructing the ICMS-evoked EMG vectors at a given site, the synergy *ν*I *n* most "evoked" at the site was the one with the largest weighting coefficient *w*I *n*, averaged over ICMS trains (**Figure 2A**, showing results for site XI in **Figure 1**, *right*). To determine which synergy was instead most "encoded" by a unit, the unit's firing rate profile was correlated against the task-related synergy scaling profiles *w*G *n* , and the largest positive correlation was identified. Results were insensitive to the use of linear (Pearson) or rank (Kendall or Spearman) correlation. To determine which task-related synergy was most evoked at a recording site, the site's ICMS-evoked EMG patterns were decomposed into combinations of synergies, which in turn were matched to taskrelated synergies (as described above). Over units, we counted the frequency at which the most-encoded synergy was the same as the most-evoked synergy at the electrode, defined as above. The chance frequency of such matches was 1/6 = 17%, as the analysis considered the 6 muscle synergies both evoked by ICMS (in the case of G2, **Figure 2A**, *left*) and observed in the task data (**Figure 2B**). Results of this χ2 test were evaluated for significance at a *p* < 0.05 threshold.

# **RESULTS**

In a recent study (Overduin et al., 2012), we examined ICMSevoked movements and muscle activity in two rhesus macaques ("G1" and "G2"). EMG data were recorded from 15–19 electrodes chronically implanted in muscles of the shoulder, arm, and hand. As we reported (Overduin et al., 2012), ICMS within MI and dorsal and ventral premotor (PMd and PMd) cortex (**Figure 1**) appeared to drive the forelimb toward an invariant, site-specific posture, and at the same time to replace voluntary muscle activity with an invariant, site-specific tonic EMG pattern (as in Griffin et al., 2011).

After aggregating ICMS-evoked EMG data over multiple stimulation sites, these pooled data could be decomposed into

### **FIGURE 2 | Microstimulation-evoked muscle activity and**

**reconstruction by synergies. (A)** The EMG activity evoked during the first 0.15 s of ICMS at site XI from **Figure 1**, shown both in gold-colored profile (*far right*) and as the integrated mean ± SE of this activity over stimulation trains (*second from right*). (The vertical 20-μV scale bars to the right of the EMG profiles give the relative level of activation in each channel.) The mean activity could be reconstructed by six ICMS-evoked synergies *ν*I *n* (*left*), with the largest weighting given to *ν*I 4. **(B)** Of eight muscle synergies *ν*G *n* derived from an object grasping behavior, synergies *ν*G 1 -*ν*G 6 could be matched in structure to ICMS-evoked synergies *ν*I 1-*ν*I 6.

combinations of a reduced set synchronous synergies using NNMF (Lee and Seung, 1999; Tresch et al., 1999). Each of the ICMS-derived synergies *ν*I *n* (*n* = 1,. . . , *N*I ) captures a pattern of synchronous firing over muscles. We found that *N*I = 6 (G2) or 7 (G1) synergies were sufficient to reconstruct ≥ 95% of the variability in the ICMS-evoked EMG data (Overduin et al., 2012). Such data reconstruction is exemplified in **Figure 2A** for G2's ICMS site XI, in MI. The figure depicts the ICMS-evoked activity at this site in the form of: the time-varying, trial-averaged EMG signal (*far right*); an ICMS-evoked EMG vector integrating this per-ICMS activity over time (*middle*); and as a vector sum of ICMS-derived synergies, each scaled by the weighting coefficient *w*I *n* (*left*). For this site, the ICMS-evoked EMG activity was dominated by synergy *ν*I 4, as shown.

We also studied a manual behavior performed by these subjects (Overduin et al., 2008, 2012), in which they had to reach for an object presented in a well, grasp it, and then carry it to the opposing well. To elicit a variety of hand postures and forces, the objects included 25 spheres, cubes, and cylinders of different dimension. Pooled over multiple days, the EMG data could be decomposed into combinations of grasping-related synergies *ν*G *n* (*n* = 1,. . . ,*N*G). We found that *N*G = 8 (G2) or 10 (G1) synergies were sufficient to reconstruct ≥ 95% of the variability in the grasping-related EMG data (Overduin et al., 2012). For each monkey, 6 of these synergies could be matched uniquely to one of the ICMS-derived synergies (**Figures 2A** vs. **2B**) using a greedy search procedure. **Figures 2B**, **3A** illustrate the synergies found for one animal and their reconstruction of average EMG activity.

Synergies could provide the animal with a mechanism to continuously and efficiently control its muscles, by specifying a task-specific amplitude coefficient time course *w*G *n* for each of the synergies (**Figure 3A** *top* vs. *bottom*). The involvement of motor cortex in this control is circumstantially suggested by the correspondence in synergy structure whether these are derived from ICMS-evoked or voluntarily-generated EMG data (**Figures 2A** vs. **2B**). Here, we sought more direct evidence for the role of motor cortex. In particular, we looked for evidence that the muscle synergies evoked by ICMS (and matched to those generated voluntarily) were also encoded by single motor cortical units near the ICMS electrodes.

A positive example of such a correspondence is shown in **Figure 3A** (*top*) and **Figure 3B** (*top*), for monkey G2. The MI unit shown in **Figure 3B** is at site XI, the location of which is highlighted in **Figure 1** (*right*), and the ICMS-evoked EMG activity of which is shown in **Figure 2A**. The single unit recorded at this site exhibited a strong burst of activity just prior to the time when the object was retrieved from its origin well, and was largely quiet during the following carry and release movements

(**Figure 3B**). (These modulations were captured in the unit's firing rate, derived from spike counts binned at 50 ms and smoothed with a 50-ms Gaussian kernel.) Of all the grasping-related synergies *ν*G 1 -*ν*G 6 , the one whose activation profile was most positively correlated with that of the unit was *ν*G 4 (**Figure 3A**, *top*), which appeared to capture the coactivation of forearm flexor muscles (**Figure 2B**). And of all the ICMS-evoked synergies *ν*I 1 − *ν*I 6, it was the matching synergy *ν*I 4 that was most dominant in reconstruction of the evoked EMG activity.

How consistently did the muscle synergy most encoded by a unit match the synergy most evoked by ICMS at the same site (as in the foregoing example)? For this analysis we simply counted the number of MI, PMd and PMv units for which the grasping-related synergy most strongly encoded by the unit matched the primary muscle synergy evoked by ICMS at the electrode. (This analysis is restricted to those ICMS sites at which units were also recorded; these are highlighted in **Figure 1**.) As shown in **Figure 4**, the actual frequencies were significantly higher than the 1/6 = 17% chance frequency [26 of 94 units, or 28%; χ2 (1) = 6.82, *p* < 0.01], as was true for G1 alone [22/83 = 27% of units; χ2 (1) = 4.82, *p* < 0.05] and supported by a trend among G2's smaller population of units [4/11 = 36% of units; χ2 (1) = 2.56, *p* = 0.07]. (Note that the "1st most encoded synergy" in this plot is not necessarily *ν*G 1 , but instead whichever synergy *ν*G *n* was most strongly correlated with the unit's firing profile, e.g., *ν*G 4 for the unit shown in **Figure 3B**. Similarly, the "most evoked synergy" is whichever synergy *ν*I *i* was dominant in reconstruction of the EMG activity evoked at the unit's electrode. In general, *ν*I 3, *ν*I 4, and *ν*I 5 appeared to be most commonly evoked by ICMS, being the synergies "evoked" at 18/83, 28/83, and 20/83 of monkey G1's units, respectively, and by 3/11, 3/11, and 4/11 of G2's units.)

While the qualitative pattern for G2 followed very closely that of G1 (**Figure 4**), it likely failed to reach significance due to insufficient sampling of units (11 units, vs. 83 for G1). Also, these trends were weakened by inclusion of PMv units. Considering *only* MI and PMd units, the fraction of cases in which the grasping-related synergy most strongly encoded by a unit matched the synergy most strongly evoked by proximal ICMS was significantly higher than 1/6 in G2 [4/10 = 40% of units; χ2 (1) = 3.27, *p* < 0.05] as well as G1 [20/73 = 27% of units; χ2 (1) = 5.04, *p* < 0.05], and in both animals combined [24/83 = 29%; χ2 (1) = 7.47, *p* < 0.01].

While the synergies *ν*I *i* extracted from ICMS-evoked EMG activity did cluster non-uniformly on the cortical surface (Overduin et al., 2012), we observed no tendency of the subset of sites at which similar synergies were both encoded and evoked to be topographically grouped on the cortical surface.

# **DISCUSSION**

Our results suggest that direct stimulation of patches of motor cortex generates synergistic combinations of muscle activity that are weighted toward synergies encoded at the stimulating electrode. A common motor cortical substrate appears to be activated both by endogenous currents during voluntary behaviors and by exogenous currents introduced by electrical microstimulation. In either case, these currents appear to modulate the amplitude, not structure, of a shared set of downstream muscle synergies. Our earlier work also suggests that each synergy (and the posture reached through its tonic activation) may be represented by a non-uniform map over the motor cortical surface (Overduin et al., 2012).

It may be contested that electrical currents injected into the CNS via microstimulation elicit non-natural patterns of neural activity (even if the downstream muscle activity is resolvable into well-organized movements and natural muscle synergies). The biological basis of ICMS, and even the extent of cortex activated, remain poorly understood (Butovas and Schwarz, 2003; Tolias et al., 2005), despite recent studies using single-cell recording, behavioral methods, functional magnetic resonance imaging (Tehovnik et al., 2006), and two-photon calcium optical imaging (Histed et al., 2009). Moreover, these investigations of ICMS have largely been limited to more conventional, short-train and/or subthreshold ICMS rather than the form used here. In applying ICMS in motor cortex, researchers have traditionally used short train durations (typically 25–70 ms) and sub- or per-threshold currents (typically 10–60μA; e.g., Asanuma and Rosén, 1972; Sato and Tanji, 1989; Donoghue et al., 1992) in order to map the overt response of cell populations including corticomotoneuronal cells. Such studies neither sought nor reported convergent movement responses of the sort found here and by others (Graziano et al., 2002a,b), and have naturally aroused some controversy (Strick, 2002).

Together with Overduin et al. (2012), our findings indicate that movements, muscle patterns, and cortical activations evoked by relatively long-train ICMS can be related to those observed in natural behavior. Indeed, it can be argued that a physiologically realistic model of motor activation *requires* ICMS with relatively long trains, as well as suprathreshold currents and intermediate pulse frequencies (optimally 80–140μA and 80–140 Hz, in the case of primate MI forelimb-area ICMS; Van Acker et al., 2013). Graziano and coworkers (Graziano et al., 2002a) have emphasized that ∼500-ms stimulation trains approximate the time scale of natural movements like primate reach and grasp (Georgopoulos et al., 1986; Reina et al., 2001). Even in spinalized amphibians, in whom convergent force patterns were first observed (Mussa-Ivaldi et al., 1990; Giszter et al., 1993; Loeb et al., 1993), these movements were evoked by similarly long-train (typically 300-ms) intraspinal microstimulation. Researchers using ICMS to study oculomotor and somatosensory systems have also used relatively long trains. For instance, 400-ms trains were applied to the arcuate sulcus to replicate the time scale of typical regular head movements (Freedman et al., 1996), and 500-ms trains were applied to primary somatosensory cortex to mimic a tactile stimulus (Romo et al., 1998). While trains shorter than ∼500 ms elicit truncated movements when delivered to primate MI (Graziano et al., 2002a; Van Acker et al., 2013), in Overduin et al. (2012) we show that even 150 ms of ICMS yields convergent movements with a predictable equilibrium point and synergistic muscle activation.

For the majority of units we sampled, firing rate profiles were more correlated with muscle synergies *other* than the mostevoked synergy (**Figure 4**, non-shaded bars). There are many ways to account for these other units. For example, units may encode other continuously-controlled quantities like the position, velocity, or acceleration of extrinsic effectors or intrinsic joints, or the forces underlying these kinematics, or the muscle contractions determining these dynamics (Carmena et al., 2003; Morrow and Miller, 2003; Schieber and Rivlis, 2007; Stark et al., 2007; Velliste et al., 2008; Ganguly and Carmena, 2009; Hochberg et al., 2012; Collinger et al., 2013). Motor cortical encoding of purely extrinsic variables may be unlikely based both on theoretical grounds (Todorov, 2000; Paninski et al., 2004) and on empirical evidence that neurons' intrinsic muscle fields are more stable than their tuning to extrinsic dimensions like hand direction (Morrow et al., 2007), with such tuning observed to fluctuate even during stable within-session behavior (Carmena et al., 2005; Rokni et al., 2007; cf. Chestek et al., 2007). However, the purpose of this report was not to determine whether motor cortical units encode muscle activity, synergy coefficients, or any other motor variable. Instead, we simply sought to test whether the units were more likely than chance to encode those synergies evoked by ICMS at the same electrode.

Another way to account for units which appeared to encode different synergies than those evoked by nearby ICMS is to recognize that the number of units affected by ICMS current far exceeds the handful sampled at an electrode. With the stimulation currents needed to elicit complex movements (≤ 100μA in this report), cortex is no doubt activated far outside the 100-μm radius assumed for ∼ 10μA stimulation (Ranck, 1981). Such suprathreshold trains are likely to recruit other cells beyond those immediately next to the stimulating electrode through synaptic connections (Histed et al., 2009), and indeed may be required to transsynaptically activate non-direct connections between motor cortex and the spinal cord (Strick, 2002).

Another possibility to be explored in future work is that populations of units may specify not only continuous variables like synergy amplitude, but also sequential recruitment of muscles (e.g., via synergies) in the form of discrete "motor programs" (Keele, 1968; Polit and Bizzi, 1978, 1979; Georgopoulos et al., 1983). Elsewhere these programs have been referred to as "timevarying synergies," in distinction from the "synchronous synergies" discussed in the present report (Kargo and Nitz, 2003), and have been extracted from muscle data using a modified form of NNMF (d'Avella et al., 2003, 2006; Overduin et al., 2008). A dynamical systems approach may also be appropriate for capturing time-varying patterns within population-level neural and muscular data (Churchland et al., 2012).

Besides studies of voluntary neural and muscular dynamics, microstimulation studies are also broadly consistent with the idea of discrete movement encoding. Intraspinal microstimulation of sufficient duration generates forces that tend to drive a limb to particular postures or through sequences of postures, in frogs (Mussa-Ivaldi et al., 1990; Giszter et al., 1993; Kargo and Giszter, 2000), cats (Lemay and Grill, 2004) and rats (Tresch and Bizzi, 1999). Sufficiently-long ICMS trains applied to mammalian motor cortex (Ward, 1938; Graziano et al., 2002a, 2004a, 2005; Stepniewska et al., 2005, 2011; Ramanathan et al., 2006) can evoke complex multijoint forelimb behaviors with multiple phases (such as reach, grasp, then retraction) and invariant endpoints. The stimulation-evoked convergent forces (Giszter et al., 1993; Kargo and Giszter, 2000), invariant endpoints (Graziano et al., 2004a; Graziano and Aflalo, 2007), and bell-shaped speed profiles (Graziano et al., 2005) all tend to overlap with motions and postures found in subjects' natural behavior. Together with the present results, these studies indicate that microstimulationevoked movements are not artifactual, and indeed can provide insights into natural movement planning.

# **ACKNOWLEDGMENTS**

This work benefited from the assistance of Margo Cantor, Charlotte Potak, Jinsook Roh, and Sylvester Szczepanowski. The project was supported by fellowships to Simon A. Overduin from the Canadian Institutes of Health Research and Dystonia Medical Research Foundation, and a National Science Foundation grant EFRI-1137267 to Jose M. Carmena.

# **REFERENCES**

Giszter, S. F., Mussa-Ivaldi, F. A., and Bizzi, E. (1993). Convergent force fields organized in the frog's spinal cord. *J. Neurosci.* 13, 467–491.

circuits underlying motor behavior. *Proc. Natl. Acad. Sci. U.S.A.* 108, E725–732. doi: 10.1073/pnas.1109925108

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 17 November 2013; accepted: 09 February 2014; published online: 05 March 2014.*

*Citation: Overduin SA, d'Avella A, Carmena JM and Bizzi E (2014) Muscle synergies evoked by microstimulation are preferentially encoded during behavior. Front. Comput. Neurosci. 8:20. doi: 10.3389/fncom.2014.00020*

*This article was submitted to the journal Frontiers in Computational Neuroscience. Copyright © 2014 Overduin, d'Avella, Carmena and Bizzi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Motor cortical regulation of sparse synergies provides a framework for the flexible control of precision walking

# *Nedialko Krouchev and Trevor Drew\**

*Groupe de Recherche sur le Système Nerveux Central, Département de Physiologie, Université de Montréal, Montréal, QC, Canada*

### *Edited by:*

*Andrea D'Avella, IRCCS Fondazione Santa Lucia, Italy*

### *Reviewed by:*

*Simon Giszter, Drexel Med School, USA Paul S. G. Stein, Washington University in St. Louis, USA*

### *\*Correspondence:*

*Trevor Drew, Groupe de Recherche sur le Système Nerveux Central, Département de Physiologie, Université de Montréal, CP 6128 Succ. Centre-ville, Montréal, QC H3C 3J7, Canada e-mail<EMAIL_ADDRESS>
We have previously described a modular organization of the locomotor step cycle in the cat in which a number of sparse synergies are activated sequentially during the swing phase of the step cycle (Krouchev et al., 2006). Here, we address how these synergies are modified during voluntary gait modifications. Data were analysed from 27 bursts of muscle activity (recorded from 18 muscles) recorded in the forelimb of the cat during locomotion. These were grouped into 10 clusters, or synergies, during unobstructed locomotion. Each synergy was comprised of only a small number of muscles bursts (sparse synergies), some of which included both proximal and distal muscles. Eight (8/10) of these synergies were active during the swing phase of locomotion. Synergies observed during the gait modifications were very similar to those observed during unobstructed locomotion. Constraining these synergies to be identical in both the lead (first forelimb to pass over the obstacle) and the trail (second limb) conditions allowed us to compare the changes in phase and magnitude of the synergies required to modify gait. In the lead condition, changes were observed particularly in those synergies responsible for transport of the limb and preparation for landing. During the trail condition, changes were particularly evident in those synergies responsible for lifting the limb from the ground at the onset of the swing phase. These changes in phase and magnitude were adapted to the size and shape of the obstacle over which the cat stepped. These results demonstrate that by modifying the phase and magnitude of a finite number of muscle synergies, each comprised of a small number of simultaneously active muscles, descending control signals could produce very specific modifications in limb trajectory during locomotion. We discuss the possibility that these changes in phase and magnitude could be produced by changes in the activity of neurones in the motor cortex.

**Keywords: locomotion, motor cortex, voluntary gait modifications, cat, synergy**

# **INTRODUCTION**

The question of modularity within the locomotor control system has a long history. From the original proposition of Graham Brown (1911, 1914) that two half-centres are responsible for generating locomotion grew a body of experimental work to determine the neuronal basis of the alternating rhythmical activity observed during locomotion (Jankowska et al., 1967a,b; Lundberg, 1981). This work took a new and, initially, controversial, direction when Grillner postulated the existence of a central pattern generator (CPG) for locomotion (Grillner and Zangger, 1975, 1979; Grillner, 1981). Subsequent studies have shown that this CPG in the spinal cord has the intrinsic ability to generate an intricate and complex pattern of locomotor activity (Pearson and Rossignol, 1991). From the original concept of the CPG as a single neuronal entity (or network) has grown the idea that the CPG may in fact be considered as a series of modules that are interconnected and provide the capacity to produce a rich behavioral repertoire, involving flexible and coordinated activity around multiple joints.

A key evolution in this respect was the formulation of the concept that the CPG is comprised of a number of unit pattern generators (Grillner, 1981; Grillner and Wallen, 1985). Grillner suggested the existence of four pairs of unit pattern generators responsible for producing the rhythmical pattern of activity in the hip, knee, ankle, and toe muscles. He proposed that relatively simple command signals to modify the connections between the hip and knee modules, for example, could easily change the pattern of muscle activity required for forward progression into that required for walking backwards.

In the turtle spinal cord, such an organization of interconnected modules has been demonstrated electrophysiologically to generate several different scratch patterns depending on the location of the offending stimulus (Berkowitz and Stein, 1994; Stein and Smith, 1997). Additional evidence for modularity in the turtle has come from the studies of deletions, in which bursts of activity in some populations of interneurones and their associated motor pools are absent in some scratch cycles while others persist (Stein and Daniels-McQueen, 2002; Stein, 2008). However, the neuronal mechanisms leading to modularity within the mammalian spinal cord are more difficult to study because of the size and complexity of the neuronal networks. As in the turtle, some evidence for modularity, although not necessarily for independent burst generators, has also come from the study of deletions (Grillner and Zangger, 1979; Jordan, 1991; Smith et al., 1995; Lafreniere-Roula and McCrea, 2005; Zhong et al., 2012). These deletions are observed in flexor and extensor motoneurones and, in many cases, have no effect on cycle timing, supporting views that rhythm generation (defining the basic rhythmicity of the step cycle) and pattern generation (defining the spatio-temporal organization of the muscle bursts within the step cycle) are separate (Lennard, 1985; Koshland and Smith, 1989). These results have led to a particularly interesting computational model of locomotion consisting of a rhythm generator and distinct pattern generators (Rybak et al., 2006; McCrea and Rybak, 2007; Zhong et al., 2012).

As an extension of these ideas, and particularly on the basis of Grillner's unit CPG model, we suggested (Drew, 1991a) that the existence of such modules could provide a substrate by which the motor cortex could exert a precise control over the magnitude, duration and relative timing of specific muscles groups while at the same time ensuring that these changes are appropriately integrated into the locomotor cycle. Similar ideas were incorporated into computer simulations designed to determine how supraspinal command signals might interact with spinal unit CPGs in a human model of locomotion (Taga, 1995, 1998).

More recently, the idea of modularity within the spinal circuits has been developed to incorporate the idea of muscle synergies. This development, which owes much to the studies of Bizzi and his collaborators (Bizzi et al., 1991; Tresch et al., 1999, 2002; d'Avella et al., 2003; d'Avella and Bizzi, 2005), posits that the nervous system produces complex movements by combining the activity of a limited number of synergies, 4–6 in most studies. In brief, synergies are defined mathematically, by using decomposition methods such as non-negative matrix factorization (NNMF), as a matrix of weights that differentially activate *all* of the muscles involved in producing a movement. By modifying the magnitude and the phase of activity of each synergy, a wide range of movement patterns can be produced (Tresch and Jarc, 2009). Such synergies have been described during locomotion, scratching and swimming in the frog (Giszter et al., 1993; Saltiel et al., 2001; Cheung et al., 2005), during postural compensation to perturbation in cats and humans (Ting and Macpherson, 2005), during human locomotion (Ivanenko et al., 2004; Lacquaniti et al., 2012) as well as in reaching movements in primates and humans (d'Avella et al., 2006, 2008; Overduin et al., 2008). Muscle synergy analysis has also been used to study the deficits in movement after stroke (Cheung et al., 2009, 2012; Clark et al., 2010). In some of these experiments, and particularly those performed in spinal animals, these synergies have been suggested to form the basis of unit burst generators of the type proposed by Grillner (Hart and Giszter, 2004; Cheung et al., 2005).

In our own studies of synergies (Krouchev et al., 2006; Drew et al., 2008a), we have taken a different approach in which synergies are defined using more classical, physiological methods. In our approach, a synergy is defined as a group of muscles that are synchronously activated such that the period of activity during locomotion begins and ends simultaneously in all muscles in a synergy. By applying a custom clustering algorithm we were able to identify 11 synergies in the forelimb of the cat during unobstructed locomotion, with nine of these occurring during the swing phase. These synergies differed from those identified most commonly by using mathematical decomposition methods (see above) in that they were both more numerous and contained only a small proportion of the total number of muscle bursts recorded. We refer to these as sparse synergies to differentiate them from those synergies comprised of all periods of muscle activity (see preceding paragraph). Our suggestion, as in our original work on this issue (Drew, 1991a), is that the motor cortex modulates these sparse synergies in order to modify limb trajectory during locomotion (Drew et al., 2008a,b).

The question remains, however, as to *how* these synergies are modified during voluntary gait modifications and whether the modifications in synergies that are observed are compatible with the discharge activity of neurones in the motor cortex recorded during the same behaviors. To answer that question, we have extended our previous analysis to examine how synergies are modified when cats step over obstacles. We consider two conditions. First we examine how synergies are modified in the lead (when a given limb is the first to step over an obstacle) and the trail (second to step over the obstacle) condition, for which our previous studies showed differential patterns of muscle activation related to the biomechanical requirements of the gait modification (Drew, 1993; Lavoie et al., 1995; McFadyen et al., 1999). Second, we examine how the synergies are modified when the cats step over obstacles of different shapes and sizes (Drew, 1988, 1991b, 1993). We then discuss the implications of these results for the cortical control of voluntary gait modifications.

# **METHODS**

The methods used in this study were either identical or similar to those previously detailed (Krouchev et al., 2006) and will only be briefly described here. All animals used in this study were originally used for other studies (Drew, 1993; Stapley and Drew, 2009; Yakovenko et al., 2011). Details concerning the methods used for animal training and implantation methods in the different animals used in the current manuscript can be found in those manuscripts. Details concerning the analytical methods are found in Krouchev et al. (2006).

# **TRAINING**

Cats were trained to walk on a treadmill at 0.35–0.5 m.s−1, first in the absence of any obstacles (unobstructed locomotion) and then in the presence of 1 or 2 obstacles attached to the treadmill belt. In these experiments, the obstacles always moved at the same speed as the treadmill. The obstacles were visible to the cat a minimum of two steps before the step over the obstacle. Most data were obtained from studies in which the cat stepped over a cylindrical obstacle of 10 cm cross-section. Additional data were obtained during steps over a high obstacle (13 cm high, 1 cm wide), a small high obstacle (7.5 cm high, 1 cm wide) and a wide obstacle (2.5 cm high, 15 cm wide) (see Drew, 1993 and **Figure 8E**).

# **SURGERY**

After the cats were trained, they were prepared for surgery under general anesthesia and in aseptic conditions following the protocols approved by the animal ethics committee at the Université de Montréal and according to the recommendations of the Canadian Council for the Protection of Animals. Details for different animals can be found in the references above. Heart rate and temperature were monitored during surgery and anesthesia level was regularly verified by testing for a corneal reflex. Solutes were administered throughout the surgery and analgesics were provided prior and subsequent to the surgery.

In all animals, pairs of Teflon-insulated, braided stainless steel wires were implanted into multiple muscles of the left forelimb. Wires were led sub-cutaneously to a 51 pin connector attached to the cranium of the cat. Although other surgical procedures were practiced on these cats, they are not reported here, as they are not relevant to the data illustrated. Details of these supplementary procedures can be found in the published manuscripts (see above). All data presented in this manuscript were obtained several weeks following the surgical procedures.

### **PROTOCOL**

Electromyographic (EMG) activity was normally recorded along with cell activity during locomotion. Data were generally recorded during unobstructed locomotion and then during periods of locomotion when cats stepped over 1 or 2 obstacles attached to a treadmill belt. EMG data were filtered at 100 Hz to either 450 or 500 Hz and amplified by a factor of 1–10 K to produce a signal of ∼1 volt. Data were digitized at 1 KHz. Video recordings of the locomotion were obtained simultaneously with the EMG data and synchronized with the aid of a SMPTE (Society of motion pictures and television engineers) digital time code.

## **ANALYSIS**

The video recordings were initially screened for sections of stable locomotion in which the cat maintained its position in the center of the treadmill. As such, we chose sections of data in which EMG activity was uniform from cycle to cycle with no evidence of intermittent changes. We then manually marked the onset and offset of each burst of EMG activity within these sections (**Figure 1A**) and classified them as either steps over the obstacle with the left leg leading or trailing, or as unobstructed locomotion (either no obstacles attached to the treadmill or at least two steps before a step over an obstacle; see Drew, 1993). A step cycle was defined as the time between two successive periods of activity in either the brachialis (Br; **Figures 1**–**7**) or the cleidobrachialis (ClB; **Figure 8**) muscle, each of which becomes active at approximately the onset of swing. The step cycle was normalized to unity (1.0) and the onsets and offsets of periods of activity in each selected muscle were then expressed as a proportion, or phase, of the step cycle. Measured events occurring after the onset of the ClB or Br were given positive values while those occurring before ClB or Br onset were given negative values. Data were plotted in phase space in which the phase of offset of a given burst was plotted as a function of its onset (**Figures 1C**, **2A**).

### *Associative clustering*

Trial data points are formed pairwise in the phase plane (*x* = *onset*, *y* = offset) and one or more bursts may be *associated* in clusters as described in Krouchev et al. (2006). It is important in this respect to realize that the data points for each muscle

**FIGURE 1 | Data selection and clustering algorithm. (A)** Untreated data during locomotion from the brachialis (Br), Cleidobrachialis (ClB), and Biceps brachii (Bic) muscles. Upward and downward directed arrows indicate onset and offset of muscle burst activity, respectively. All activity was synchronized to the Br and burst onset activity is defined as phase = 0.0. A step cycle is defined as the time between two successive bursts (phase = 1.0). Activity in other bursts is defined with respect to Br onset. The red arrows indicate data illustrated in **(B,C)**. The green arrows define the first burst of activity in Bic [Bic(1) in **Figure 2**]. **(B)** The rectangles define 1SD of the mean onset and offset of the activity of each muscle (P1–P3, individual data illustrated in **C)**. A vector is drawn from the center of each rectangle to the vertex closest to the center of the nearest neighbor (thick lines, Sn). The vector distance (dotted line: Dn) is then calculated as described in the text. **(C)** Bivariate ellipses are drawn around the centroid of each cluster (see text). Same dataset as in **Figure 2**.

are defined as belonging to that muscle and cannot be divided among different clusters. In the clustering algorithm, it is assumed that each burst is fully described by its centroid—i.e. the mean (onset, offset) phase vector, and the associated pair of standard deviations (SD's: σ*X*, σ*Y*; assumed uncorrelated). Hence, each burst is assumed equivalent to the encompassing rectangle *[*−σ*X*, σ*X]* × *[*−σ*Y*, σ*Y]* around the centroid mean values *[X*, *Y]* (**Figure 1B**).

Sufficiently overlapping rectangles form clusters, or synergies (Krouchev et al., 2006). Thus the algorithm associates sets of data points from different muscles rather than dissociating them as do most clustering methods.

First, the individual data points are tested for outliers using Rosner's test (Rosner, 1983) with a 2 SD margin. For each burst *i*, the phase-plane vectors of mean phase:

$$P\_i = (X\_i, Y\_i) = mean \left( onset\_i, offset\_i \right)$$

and of standard deviations

$$S\_i = (\sigma X\_i, \sigma Y\_i) = std\left(onset\_i, offset\_i\right).$$

are calculated. The distance between two bursts *i* and *j*, is expressed as the phase-plane vector (**Figure 1B**):

$$D\_{\vec{\eta}} = P\_i - P\_j$$

Two bursts are considered part of the same cluster when:

$$q\_{i\bar{j}} = \max\left\{ \left( P\_{\bar{i}} \pm \mathcal{S}\_{\bar{i}} \right) D\_{i\bar{j}} + \left( P\_{\bar{j}} \pm \mathcal{S}\_{\bar{j}} \right) D\_{i\bar{j}} - ||D\_{i\bar{j}}||^2 \right\} / ||D\_{i\bar{j}}|| > 0 \tag{1}$$

where the maximum is taken over all possible combinations (**Figure 1B**), the products are the dot (scalar) vector products, and || · || is the usual Euclidean norm in the xOy phase plane.

For all possible pairs of bursts *i* and *j* (*i* = 1, 2 ..., *n*, *j* = *i* + 1,..., *n*)— where *n* is the total number of bursts to classify, Equation (1) is verified and thence a boolean *n* × *n* uppertriangular adjacency matrix **B** is obtained (**B***ij* = 1 when **q***ij* > **0**, and zero otherwise). The latter adjacency matrix is then the input to a Matlab implementation of an algorithm due to Press et al. (1992)to determine the equivalence classes—i.e., the clusters. The final cluster numbers are rearranged so that their mean onsets proceed in an ascending order.

The clustering proceeds in two sub-stages. First, we cluster the bursts coming from the first (and most elaborate and ample) subset of the data—coming from the same animal. In the second sub-stage, additional bursts (coming from additional animal data sets) are either allowed to join an existing cluster, or to form a new one. In cases in which the existing cluster has more than one member, the new candidate burst will join if it overlaps with at least two bursts from the existing cluster. It should be noted that this associative clustering method is robust to even relatively large changes in the centroids of EMGs making up a cluster. In a previous publication (Krouchev et al., 2006; supplemental information), we tested the stability of the methods by introducing random jitter to the centroids of the muscles making up a cluster. In 970/1000 cases in which each centroid was displaced by <0.6 SD, there was no change in cluster composition compared to that obtained using the actual data.

### **CONSTRAINED vs. UNCONSTRAINED CLUSTERS**

In the analysis in our original papers (Krouchev et al., 2006; Drew et al., 2008a) the algorithm was always unconstrained in that muscles were formed into clusters, and therefore synergies, simply on the basis of the rules summarized in the preceding paragraph. A similar approach was used to obtain the synergies active during unobstructed locomotion (**Figure 2**) as well as during the lead and trail conditions of the voluntary gait modifications. However, because of the addition of periods of muscle activity that occurred only during the gait modifications, and because of the relative changes in phase of some muscles, there were small changes in the composition of some of the synergies. This makes a direct comparison of changes in the phase and magnitude of synergies during all three conditions problematic. We, therefore, also used a constrained analysis in which muscles bursts were confined to the same clusters identified during unobstructed locomotion. This allowed us to directly compare the changes in phase and magnitude of synergies comprised of the same muscles in all three conditions. Further details of the approach and of its limitations are provided in the Results and Discussion.

## *Direct component analysis (DCA)*

For each cluster we derive a temporal activation profile, which is labeled a *direct* component (DC). The latter contrasts to the muscle activation components in the literature, which most often are obtained through abstract mathematical decomposition.

DCs are derived directly from the corresponding overall **(**onset, offset**)** phase statistics, which describe the periods of activity of the muscles (or bursts) forming the cluster. For each cluster, we calculate the 2 *marginal univariate* Gaussian probability density functions (pdf) for the onset and the offset phases. For each cluster *k*, the onset/offset pdf is, respectively, *N*(*Xk*, σ*Xk*) and *N*(*Yk*, σ*Yk*) and activation of the muscles forming this cluster is assumed to span the phase interval (*Xk* − 3σ*Xk*, *Yk* + 3σ*Yk*).

We further assume that the overall shape of the temporal activation profile for cluster *k* is captured by Gaussian basisfunctions. Hence its DC *uk(t)* is described by:

$$\mu\_k(t) = \begin{cases} \exp\left(-\left(t - Z\_k\right)^2 / 2\sigma\_1\right), & t < Z\_k\\ \exp\left(-\left(t - Z\_k\right)^2 / 2\sigma\_2\right), & t > Z\_k \end{cases}$$

where

$$\begin{aligned} Z\_k &= (X\_k + Y\_k)/2, \sigma\_1 = (Z\_k - X\_k)/3 + \sigma X\_k, \text{ and} \\ \sigma\_2 &= (Y\_k - Z\_k)/3 + \sigma Y\_k \end{aligned}$$

### *Cluster statistics and representation*

In this paper we use a correlated bivariate-normal model. The data scatter for burst *k* is captured by the pdf:

$$f\left(\mathbf{x},\mathbf{y}\right) \sim \exp\left(-z/2\sigma\right) \tag{2}$$

where:

$$\begin{array}{rcl} z = \mu^2 + \nu^2 - 2 \ \mu\nu\rho & \nu = 1 - \rho^2\\ \mu = (\varkappa - X\_k) \ / \sigma X\_k & \nu = \left(\nu - Y\_k\right) / \sigma Y\_k \end{array}$$

Hence, each cluster is represented by an ellipse of the form:

$$
\mu^2 + \nu^2 - 2\,\mu\nu\rho = 1 - \rho^2 \tag{3}
$$

The ellipses' skewing and orientation depend on the value of **ρ**. Krouchev et al. (2006) assumed no (x,y) correlation—i.e., **ρ** = 0. This yielded ellipses with major axes, respectively (σ*X*, σ*Y*), parallel to the coordinate axes. For the general case of non-zero **ρ**, it may be shown that the ellipse described by Equation (3) touches the ±1*SD* encompassing rectangle exactly at the points (±**1**, ±**ρ**) and (±**ρ**, ±**1**) (**Figure 1C**).

### *Statistics*

To determine whether the phase and/or the magnitude of the synergies were significantly different in each condition we used ANOVAs across the conditions (unobstructed, lead and trail). Individual One-Way ANOVA's were performed for each of the synergies using the data for the individual periods of activity of each muscle within a synergy. Tests for difference in the phase of activity were made using the onset ∼ *N*(*Xk*, σ*Xk*), the offset ∼*N*(*Yk*, σ*Yk*) and the peak phase (*Xk* + *Yk*)/2 of activity. The nullhypothesis for each ANOVA is that condition has no significant effect on the range of phase-values in the sample. For synergies that showed a significant (*p* < 0.05) effect of the locomotion condition, we performed pair-wise *t*-tests between conditions. Individual One-Way ANOVA's were also performed to determine if the magnitude of the synergies also varied as a function of condition (see Results).

# **RESULTS**

### **DATABASE**

The database for the principal analysis in this paper is based around cat RS26, as in a previous publication (Krouchev et al., 2006) examining synergies during unobstructed locomotion. Supplementary muscles in that previous publication were obtained from 3 other cats, PCM3, MC8 and RS23. Ideally, we would have used the same database in the current manuscript. Unfortunately, however, data for voluntary gait modifications were not available from all animals. We have, therefore, complemented the data from RS26 with data from cats MC29 and RS23 for the current analysis to produce a full dataset during unobstructed locomotion consisting of 27 bursts of EMG activity from 18 muscles of the forelimb, all recorded on the same side. The 10 clusters defining the muscle synergies in this dataset during unobstructed control locomotion are illustrated in **Figure 2** and the average activity of some of these muscles can be observed in **Figure 3** (black lines). Note that the clusters defined in the present analysis are very similar to those published previously (Figure 6 in Krouchev et al., 2006). The major difference is that the current dataset produced 10 clusters of activity rather than the 11 that we defined previously. This was primarily caused by the muscles in cluster #2 forming a single cluster (including EDC, LtD and TrM) rather than being divided into 2 clusters as in our previous analysis. One other difference is the presence of a very early burst of activity in the Bic, preceding activity in all other muscles. Apart from these minor differences with our previously described dataset (Krouchev et al., 2006), the essential features of these synergies are maintained. These include: the sequential pattern of activation of the synergies; the sparse nature of each synergy; the fact that a given muscle may be represented in more than one synergy and the fact that a given synergy can include muscles acting proximally together with others acting distally (e.g., EDC and TrM in cluster #2).

# **CHANGES IN PHASE OF SYNERGIES DURING GAIT MODIFICATIONS (LEAD AND TRAIL CONDITIONS)**

During the gait modifications, there were substantial changes in the magnitude and phase of some of the bursts of EMG activity as we have also described elsewhere (Drew, 1993; Lavoie et al., 1995; Drew et al., 2008a). The changes observed in selected muscles of the current dataset are illustrated in **Figure 3** for the lead and trail condition. The EMG data for all muscles are normalized in time to the average cycle duration for the original dataset in which they were recorded, and normalized in amplitude to the largest magnitude recorded from any one given muscle. As such, changes in amplitude and relative changes in phase are clearly visible in the presentation. Several points need to be emphasized. First, in both conditions, there are changes in amplitude and phase in a number of muscles. During the lead condition, these changes are primarily expressed as changes in amplitude (e.g., Br and BrR) with fewer changes in phase, although see phase delay in EDC. In contrast, during the trail condition, there are major changes in both phase and amplitude in a number of muscles. This is particularly clear in the shoulder retractor muscles (LtD and TrM) and in the muscles acting around the wrist and digits (PrT, EDC). Last, one other major difference from the control activity is the presence of periods of activity during the gait modifications (**Figure 3**, arrows) that were not present in the unobstructed condition. The clearest examples are the supplementary periods of activity in the EDC, SpD and the Tri at the onset of swing in the trail condition. There are also small changes in some of the extensor muscles. These occur after the step over the obstacle in the lead condition and before the step over the obstacle in the trail phase.

In general, most of the changes in phase of the muscles active during swing occur during the trail condition, before the onset of swing (defined here as Br onset) while the more modest phase changes in lead occur after swing onset. As discussed previously (Drew, 1993; Lavoie et al., 1995; McFadyen et al., 1999), this difference is related to the constraints of the task. During the lead condition, swing begins with the obstacle well forward of the paw and the major requirement is to lift the limb above and over the advancing obstacle. During the trail condition, swing begins with the obstacle close to the trail paw and the major requirement is to retract the limb sufficiently to ensure that the limb is lifted away and above the obstacle as it continues to advance (see Drew, 1993).

To determine the effects of the changes in the pattern of muscle activity on synergy composition, we performed the cluster

**FIGURE 2 | Synergies during unobstructed locomotion. (A)** The results of the cluster analysis are illustrated for a total of 27 periods of EMG activity recorded from 18 muscles in 3 cats (30–70 values for each muscle burst). 10 clusters, each corresponding to a muscle synergy, are identified by color and the EMG bursts in each cluster are identified by symbol and color as shown in the key to the right. The same order of muscles is used in the key of all figures except when a change in synergy composition makes it impossible. The spatial location of each cluster is illustrated by the ellipses (thick lines). EMG onsets and offsets are referenced to the onset of activity in the brachialis (Br) muscle (phase = 0.0). Negative phase onsets indicate muscle bursts active prior to the onset of activity in the Br. **(B)** the muscles used in

this study are illustrated on stick figures of the cat's forelimb taken from Crouch (1969). Clusters are organized according to the biomechanical function and the colors of the muscles correspond to the colors of the synergies. Note that colors accorded to some muscles are different from those used in Krouchev et al., 2006 because of a slightly different composition of the synergies (see text). Muscle Abbreviations: AcD, acromiodeltoideus; BrR, brachioradialis; ClT, cleidotrapezius; ECR, extensor carpi radialis; EDC, extensor digitorum communis; ECU, extensor carpi ulnaris; LtD, latissimus dorsi; LvS, levator scapularis; PaL, palmaris longus; PrT, pronator teres; SpD, spinodeltoideus; SSp, supraspinatus; Tri, triceps bacchii, long head; TriL, triceps brachii, lateral head; TrM, teres major.

muscles during the lead (left column, red traces) and the trail (right column, green traces) condition. Average activity during gait modification is superimposed onto the activity in the unobstructed condition (black traces). Data are synchronized to the onset of the Br and the duration

analysis on the datasets for the voluntary gait modifications (lead and trail conditions) using all available bursts of EMG activity including those that were present only during the voluntary gait modifications (**Figure 4**). As explained in the Methods, we used two complementary approaches to examine changes in synergy. In the first instance, we applied the clustering algorithm to the entire dataset in the lead and the trail condition in exactly the same way as we did for the unobstructed data. We refer to this as the unconstrained condition (**Figures 4A,B**). In this approach, we are effectively asking whether the same synergies are found in all three conditions studied, unobstructed, lead and trail. In the second approach, and the one that we used for all further analyses, we constrained the synergy composition during the lead and trail conditions to be identical to that observed in the unobstructed condition (**Figures 4C,D**). In this instance, we are determining how the phase of the synergies identified during unobstructed locomotion needs to be modified in order to obtain the patterns of activity observed during the gait modifications. In this second approach, we used the adjacency matrix to determine to which cluster the bursts activated only during the gait modifications were associated. This, constrained, approach has the advantage of allowing us to provide a direct comparison of the changes in phase and magnitude that occurred in each synergy during these gait modifications.

taken from cat MC29, those without an asterisk are from cat RS26. Arrows indicate bursts of activity that were not present in the

unobstructed condition.

The results obtained by using the unconstrained algorithm show that the muscles bursts observed during the voluntary gait modifications are organized into synergies that are fully consistent with those observed in the unobstructed condition. Specifically, the analysis results in the presence of a similar number of synergies as in the unobstructed condition, with the same sequential organization as in unobstructed locomotion and with each synergy being comprised of a small number of muscles (**Figures 4A,B**).

clusters formed by applying the analysis to all bursts recorded in each condition and allowing the analysis to define the resulting clusters (unconstrained). **(C,D)** analysis performed using all bursts but with the clusters constrained to have the identical composition as in the

condition, the colors describing some clusters are different from those used in **Figure 2** because of the modification of synergies detailed in the text. In the constrained condition, the colors are identical to those used in **Figure 2**.

However, there were some differences in cluster composition, compared to the unobstructed condition as can be seen by examining the legend identifying the clusters in **Figures 4A,B**. These differences are illustrated in more detail in **Figure 5** where the synergies identified during the swing phase in the unconstrained condition (thin ellipses) are compared to those obtained in the constrained condition (thick ellipses) for the lead (**Figure 5A** and trail (**Figure 5B**) conditions (see also **Table 1**).

In the lead condition, there were only minor changes between constrained and unconstrained synergies as can be seen by the overlap of the ellipses representing the synergies obtained in the two conditions (**Figure 5A**) and by inspection of **Table 1A**. Indeed, 7 of the 10 clusters showed no changes at all in the two conditions (**Table 1A**). Among those clusters that were modified in the two conditions, the largest change was the division of cluster #2 (as defined in the unobstructed condition) into 2 clusters in the unconstrained condition. This was because of a slight phase advance of the EDC(1) burst with respect to the other three muscles in the cluster. Note, however, that even for this relatively large modification in phase, the change in centroid position of the EDC (0.059) was substantially less than the distance between the centroids of clusters #1 and #2 (0.13) and between clusters #2 and #3 (0.19) in the constrained condition. In addition, the Br was included in a cluster with the PrT and the SpD instead

illustrated as thick lines and those from the unconstrained condition as thin lines. Each symbol represents the centroid of a burst of activity in a given muscle as indicated by the key for the constrained condition (to the left of each illustration). It should be emphasized that the location of the centroid of each muscle is identical for the constrained

constrained and unconstrained clusters ({). Cluster #1 in the trail condition is not displayed because of the scale used. Note that the colors identifying each muscle and each ellipse are the same for the constrained and unconstrained conditions and identical to those used in **Figure 2**. As a result, muscles with the same colors sometimes belong to different clusters in the unconstrained condition (see key).

of being separate as in the unobstructed condition. This cluster in the unconstrained condition was only displaced by 0.008 from the location in the unconstrained condition (**Table 1A**). Despite these small changes in the unconstrained condition, it is evident that the basic elements of the synergies identified in the unobstructed condition are equally visible in the step over the obstacle.

Similar qualitative changes were seen in the trail condition as can be seen in **Figure 5B**. Again, the muscles comprising cluster **Table 1 | Phase differences between the centroids of constrained and unconstrained clusters.**

*The tables show the differences in the phases of the centroids of the constrained and unconstrained clusters as illustrated in Figures 4, 5.*

#2 in the unobstructed condition were divided, this time into 2 groups of 2. Both clusters, however, remain well separate from the next cluster in the sequence. For example, the differences of the two unconstrained clusters (#2 and #3) from constrained cluster #2 was 0.057 and 0.046, respectively (**Table 1B**), while the distance between constrained clusters #2 and #3 was 0.17, three times as large. (The distance between cluster #1 and cluster #2 was 0.29). There were also changes in the other clusters, but in all cases, these involved the transfer of one muscle from one cluster into an adjacent cluster without change in the overall sequential nature of the activation of both individual muscles and clusters. In other words, in both the unobstructed condition and during the step over the obstacle, activation of e.g., LvS, follows activation of ClB and ClT, which follow activation of Br, which follows activation of PrT and SpD . . . etc.

Overall, the underlying principle of sequential activation of sparse synergies is well supported by the analysis (see Discussion). In particular, the analysis serves to illustrate that the synergies observed during the voluntary gait modifications are fully consistent with those observed during the unobstructed condition. More specifically, they serve to illustrate that the changes in the pattern of EMG activity that results in the modified limb trajectories required to step over the obstacles can be produced by changing the phase of activity of the synergies that underlie the EMG pattern during unobstructed locomotion.

To directly compare the phases of activity of the synergies in the lead and trail conditions with those identified during the unobstructed condition, we transformed the clusters obtained by using the constrained analysis into DCs, centered on the centroids of each cluster (**Figures 6A–C**) as described in the Methods (see also Krouchev et al., 2006; Drew et al., 2008a). These DCs, for each of the three conditions, are shown superimposed in **Figure 6D**, which allows the changes in phase of activation of each group of synergies to be synthesized in an economical manner. Changes in amplitude are ignored in this display but are addressed below. During the lead condition (red traces) the changes in phase are relatively minor apart from the phase delay in clusters 7 and 8, which are responsible for the wrist dorsiflexion and the preparation of the limb for landing. There is an increase in the duration of clusters #5 and #6, which include the Bic and the ClB (see **Figure 2**), and which are responsible for transporting the limb over the obstacle. In the trail condition, the phase advances observed in the activity of individual muscle bursts (**Figure 3**) have a substantial effect on the phase activity of a number of different synergies, and particularly those active prior to the onset of swing (i.e., Clusters #1–3). Note that treating cluster #2 as separate clusters as suggested by **Figures 4A,B**, **5** does not change the basic conclusion of a phase advance of muscle activity in this phase of the step cycle.

A One-Way Anova for phase of onset, phase of offset and peak phase (three different tests) as a function of locomotion condition showed significant differences for all combinations except 1 (phase of offset for synergy #9). Pair-wise comparisons of the

three conditions likewise showed significant changes of at least one of the measured variables (onset, offset, peak) for all three comparisons (control-lead, control-trail, and lead-trail). In the control-lead comparisons, phase of onset was unchanged for four

**(A)**, trail **(B)** and unobstructed **(C)** conditions. **(D)** direct components showing the changes in phase of the synergies in the lead (red

> synergies while phase of offset was changed for all 10 synergies, reflecting the fact that most changes occur after swing onset in this condition. In the control-trail comparison, phase of offset was unchanged in three synergies. Overall these results emphasize

as suggested by the analysis in the unconstrained condition

(**Figures 4**, **5**).

that gait modifications require changes in the phase of activation of all synergies.

# **CHANGES IN AMPLITUDE OF SYNERGIES DURING GAIT MODIFICATIONS (LEAD AND TRAIL CONDITIONS)**

We measured the magnitude of the EMG activity of each burst of activity used in the analysis for the unobstructed, lead and trail conditions by integrating each 1 ms bin of activity between the measured onset and offset of the period of activity. These values were then divided by the duration of the burst (in ms) to give a value that represents the level of activity of the muscle burst in each condition, independent of any changes in duration of the muscle activity.

The results from this analysis are shown in **Figure 7A**, with each measured EMG burst being classified according to the synergy to which it is assigned. As expected from the data illustrated in **Figure 3**, muscles such as the Br, ClB and ECR(2) show a large increase in activity (up to 200%) during the lead condition while others, such as the PrT(1), and the TrM(1) show their largest changes (up to 400%) during the trail condition. In most cases, these changes were significantly different (asterisks) from control levels in both the control and lead conditions. Furthermore, in many cases, muscle amplitudes were significantly different between lead and control. Note that changes in the magnitude of activity cannot be displayed for muscle bursts that were inactive in the unobstructed situation.

In general, changes in magnitude of muscle bursts (relative to unobstructed locomotion) for clusters comprised of several muscles (e.g., clusters 2, 5, 8, and 10) were of the same sign and of similar magnitude. For example, in cluster #2 all of the muscles in the cluster showed a significant increase in activity during the trail condition. In clusters #3 and #5, each muscle shows a significant increase in the lead condition while in synergy #10, the magnitude of muscles in the lead condition is either unchanged or shows only a small increase. Only in cluster #8, does one muscle show a significant decrease in activity during the lead condition while the other five muscles either show an increase or are unchanged. A One-Way ANOVA for the changes in amplitude as a function of condition (control, lead, and trail) showed significant changes for all 10 synergies. Comparison between the pairs of conditions showed significant changes in at least seven synergies in all three conditions.

These results indicate that the voluntary gait modifications require changes in the magnitude of muscles in most synergies and that these changes in magnitude are different for the lead and trail conditions.

To provide an illustration of how control signals would need to change in both magnitude and phase in order to produce gait modifications, we used the average value for the EMG bursts in each cluster to scale the amplitude of the DCs that we illustrated in **Figure 6**. The results from this procedure are illustrated in **Figure 7B** which serves to illustrate that the gait modifications require coordinated changes in both amplitude and phase of the synergies. In the lead condition, the changes in phase or amplitude in the first two synergies (mainly responsible for paw lift) were small. Subsequently, there were changes primarily in amplitude in synergies 3–6 that are responsible for flexing the limb above the obstacle and transporting it forwards. Subsequently in synergies 7 and 8, responsible for wrist dorsiflexion and then placement of the paw on the substrate, there are increases in amplitude and a phase delay of synergy 7. Minimal changes are observed in the subsequent period of extensor muscle activity. In contrast, in the trail condition, there are both phase advances in the first 3 clusters and an increase in magnitude in clusters #2 and 3. There are minimal changes in phase or magnitude in clusters 4–6 although the slight decrease in duration of cluster #5 leads to a phase advance of the activity in clusters 7 and 8.

# **EFFECT OF OBSTACLE DIMENSIONS**

Changing the shape and size of an obstacle changes the limb trajectory as the cat steps over it (Drew, 1988). This implies that the relative magnitude, duration and timing of bursts of EMG activity in each synergy have to be precisely modulated to produce the appropriate limb trajectory. This is illustrated in **Figure 8** with a different cat (MC8) for which data were available from a more limited number of muscles but for a variety of obstacles (**Figure 8E**).

**Figure 8A** illustrates EMG activity from representative muscles during control locomotion (black trace) and during steps over a thin high obstacle (**Figure 8Eiii**) in the lead (red trace) and trail condition (green trace). As in the example illustrated in **Figure 3**, there were characteristic changes in activity with phase advance of the TrM burst prior to swing onset in the trail condition, together with an increase in EMG amplitude of the ClB and phase delay of the burst of the EDC in the lead condition. Performing the cluster analysis on the EMG activity during unobstructed locomotion (**Figure 8B**) provides the same sequential pattern of activity that was obtained with our full database, despite the smaller number of muscles bursts that were available in this cat. However, because of the smaller number of EMG bursts in this database we identified only seven synergies.

The changes in phase and amplitude in the DCs constructed from these clusters are illustrated in **Figure 8C** for the lead condition and in **Figure 8D** for the trail condition. For all of the obstacles for which data were analysed these changes are similar in form to those obtained with our main database synthesized in **Figure 7B**. Anovas showed a significant effect of both condition and obstacle on the phase and magnitude for many of the synergy. In the lead condition, the effect of obstacle was significant for synergies #3, 6, and 7. In the trail condition, there was an effect of obstacle on all seven synergies.

In general, the change in phase and magnitude were smallest for the small high obstacle (green) and largest for the very high (mauve) and round (orange) obstacles. In cluster #3, in the lead condition, for example, there is a large increase in magnitude for the very high obstacle, presumably because of the need for shoulder retraction to raise the limb above the obstacle. In clusters #4 and 5, there are large changes in amplitude and duration for the 3 largest obstacles and relatively smaller ones for the small high obstacle. In cluster #6 there are large changes for the very high and round obstacles that require the paw to be raised above the obstacles and correspondingly large phase delays in cluster #7, especially for the wide obstacle. In the trail condition, changes in amplitude of cluster #3 are particular clear for the very high obstacle and to a lesser extent in the small high and the round

### **FIGURE 7 | Change in magnitude of the synergies. (A)** the mean magnitude (+SD) of the activity in each muscle during lead (red bars) and trail (green bars) is expressed as a percentage of the activity in the unobstructed condition (100%, thin horizontal line). The magnitude is normalized to the duration of the burst (see text). Asterisks beside the bars indicate muscle bursts that were significantly different from control (*p* < 0.05). Horizontal

lines with an asterisk indicate muscle bursts that were significantly different between the lead and trail conditions. **(B)** the average change in the level of activity from each of the synergies is used to modify the magnitude of the direct components, which now provide a representation of both the change in phase and the change in magnitude of the synergies that is required to step over the obstacles in the lead and trail condition.

**of different size and shape. (A)** averaged activity of 5 representative muscles during the lead (red), trail (green) and unobstructed (black) conditions. Data are synchronized to the onset of activity in the ClB. **(B)** Cluster analysis for the database of 14 muscle bursts available from this cat. **(C,D)** direct component

each synergy when the cat steps over obstacles of different sizes and shapes in the lead **(C)** and trail **(D)** conditions. The insets amplify the traces representing synergy #6 in the lead condition and synergy #3 in the trail condition. The color code for the obstacles is represented by the color of the cats in **(E)**.

obstacle. The changes in amplitude are least for the wide obstacle (which was also the least high). These changes are needed to retract the limb above the obstacle. The changes in duration are minimal except for the ECR (cluster #6). This latter is presumably because of the greater excursion required to bring the wrist into dorsiflexion in the trail condition (Lavoie et al., 1995).

# **DISCUSSION**

We have argued in previous publications (Krouchev et al., 2006; Drew et al., 2008a) that small modifications in the activity of the sparse synergies defined by our analysis could provide a flexible manner of producing modifications in limb trajectory during walking. The results obtained in the current study support this premise by clearly demonstrating that small changes in the phase and magnitude of the synergies identified during unobstructed locomotion can produce the changes in limb trajectory required to step over obstacles of different sizes and shapes, both in the lead and the trail condition.

# **SYNERGY DEFINITION AND COMPOSITION**

Our definition of a synergy is a group of muscles that become active simultaneously and remain active for the same period of time. One of the results of this definition is that our synergies include only a limited number of muscles (hence sparse synergies) compared with those defined by mathematical decomposition that generally contain all muscles contained in the dataset (see below). Indeed, some of our synergies contain only a few, or in some cases, only one muscle (e.g., Cluster #7 in **Figure 2**). As we have argued previously (Krouchev et al., 2006), this is simply the result of recording a limited number of muscles and one would expect other wrist and digit extensors to be included in this synergy, if recorded. A correlate of using our definition of a synergy is that we define a relatively large number of synergies (10–11) compared to the 4–6 that are normally identified by the more commonly used decomposition methods (d'Avella et al., 2003; Tresch and Jarc, 2009). It is, however, interesting to note that synergies with similar compositions to those obtained using our methods can be obtained if the number of synergies defined by the linear decomposition methods is increased to 10 or 11 (Krouchev et al., 2006). Similarly, synergies comprised of a limited number of muscles can be obtained by other methods by including only those that show significant differences from lower amplitude weights (Hart and Giszter, 2004).

In this study we have taken the approach of defining our synergies on the basis of the muscle activity patterns generated during a single behavior, unobstructed locomotion, and then defining how these synergies must be modified to produce other behaviors. This differs from the more common approach in which synergies are defined on the basis of all behaviors under study (see e.g., d'Avella and Bizzi, 2005. In other words, most studies ask what parameters (weighting and phase) applied to all muscles in the dataset will result in all behaviors under study. In contrast, the question that we are asking is, how does one need to modify the magnitude and phase of the synergies obtained during unobstructed locomotion in order to obtain the patterns observed during voluntary gait modifications?

One of the results of defining synergies on the basis of a single behavior is the possibility that in addition to magnitude and phase, synergy composition may also be modified in other behaviors (see e.g., Ivanenko et al., 2005; Kargo and Giszter, 2000, see below). Indeed, we observed some changes in synergy composition in the unconstrained condition of our analysis (**Figures 4A,B**, **5**). This raises the question of whether the voluntary gait modifications are really the result of modifying synergies that define unobstructed locomotion or whether there is a real need to modify synergy composition as suggested by **Figures 4A,B**. In some cases, such as the division of cluster #2 in the unobstructed condition into two clusters (#3 and #4) in the voluntary gait modifications, there is reason to think that this might indeed reflect a true need to modify the synergies. In most other cases, however, such as the merging of clusters 4 and 5, it is probable that this reflects the variability in the measurements of the onset and offset of EMG activity together with the modifications of the adjacency table because of the addition of bursts of activity not present in the unobstructed condition. Certainly in the absence of additional independent methods to define synergies this must remain an open question. However, what is clear is that the basic sequential activation of the muscles bursts and resultant synergies observed in the unobstructed condition is maintained in the voluntary gait modifications. For example, all 4 of the muscles comprising cluster #2 in the unobstructed condition are activated subsequent to the activation of the Bic(1) and prior to the activation of the Prt and SpD. Moreover, all 4 of these muscles show similar changes in magnitude during both the lead and the trail condition (**Figure 7A**). Lastly, as illustrated in **Figure 6D**, treating these muscles as 2 clusters does not alter the basic facts that the muscles active at this time are phase advanced compared to the unobstructed condition. Similar considerations hold for the other synergies. For example, whether the Br forms part of a synergy with the PrT (**Figure 5A**) or is in a separate synergy (**Figures 2**, **5B**) does not change the fact that activity in both of these muscles occurs after the period of activity in the EDC and the TrM and before the period of activity in the ECR. Thus the basic concept of a sequential activation of a (relatively) large number of sparse synergies is fully supported by the data while the exact composition of each synergy must remain open to further investigation. Nonetheless, by maintaining the same number of synergies and the same composition, this approach has the advantage of allowing us to directly address some of the properties of the descending signals that would be required to modify gait by modifying the phase and magnitude of muscle synergies.

### **CORTICAL CONTROL OF SYNERGIES**

Our results demonstrate the existence of a relatively large number of synergies, each of which is comprised of a small number of muscles (sparse synergies) and each of which is active during only a small part of the step cycle. These synergies are activated sequentially during the step cycle and show well defined modifications during the gait modifications.

In the lead condition, the major change is the increased duration and the phase delay of the synergies active at or after swing onset. These changes lead to the increased swing duration and the delayed onset of the activity of the muscles related to paw placement. There is also an increase in the magnitude of the synergy that leads to increased elevation over the obstacles. These changes in phase and magnitude are scaled to the size and the shape of the obstacle. Wider obstacles lead to progressively increased durations and longer phase delays. Higher obstacles lead to smaller changes in duration and smaller phase delays but increased amplitude in synergies related to limb flexion. In contrast, during the trail condition, the major changes are the phase-advance and increased amplitude of the synergies active before swing onset, particularly those active at paw lift and responsible for lifting the paw above the obstacle. The data therefore illustrate how differentially modifying the activity of different synergies can lead to adaptive changes in limb trajectory that allow for avoidance of obstacles covering a relatively wide range of sizes and shapes. There seems little doubt that the phase and magnitude of these synergies could be modified to produce an infinite range of limb trajectories that could be adapted to any specific obstacle that was encountered during locomotion.

In the most parsimonious reasoning, one might expect that the activity of these synergies would be modified by control signals with similar characteristics. This is exactly what is seen in recordings from the cat motor cortex during locomotion (Drew, 1993; Drew et al., 1996, 2008a). During voluntary gait modifications, a large population of motor cortical cells (86% in Drew, 1993) show a modification in activity (increase or decrease) consisting of changes in magnitude, duration and/or the relative timing (phase of onset) of the discharge (Drew, 1993; Drew et al., 1996). Sub-populations of motor cortical cells recorded during locomotion, and particularly during gait modifications, are each activated during only a small part of the step cycle in a sequential pattern (Lavoie and Drew, 2002; Drew et al., 2008b). Examples of these activity patterns during the lead condition taken from the database used in a previous publication (Drew, 1993) are summarized in **Figure 9A** which serves to illustrate both the discrete, phasic nature of the activity patterns observed in different pyramidal tract neurones (PTNs) in the motor cortex and the sequential nature of this activation pattern. PTN1, for example, is active at the onset of the swing phase as the paw is lifted from the support surface at the onset of swing. This discharge occurs at the same time as the activity observed in muscles such as the TrM and the LtD. Subsequently, other PTNs become active coincident with the elbow flexion, characterized by activity in Br (PTN2); then with the wrist dorsiflexion, characterized by activity in the ECR (PTN3) and finally during paw placement, characterized by activity in muscles such as the EDC, TrM and the LtD (PTN4).

Particularly striking in our database (Drew, 1993) were cells whose discharge activity covaried with the first period of activity in the TrM (e.g., PTN1 in **Figure 9A**). As illustrated in **Figure 9B** these cells generally show a phase advance together with a major increase in the magnitude of their discharge in the trail condition, in the same manner as does the activity in synergy #3 (**Figures 6**, **7**), containing the TrM. Importantly, discharge activity in these cells is scaled for different obstacles in the same way as the muscle synergy.

Overall, the data from our previous recording studies suggest that there are multiple sub-populations of PTNs, each of which regulates the activity of a specific sparse synergy. We suggest that there should be as many sub-populations of cortical neurones as there are synergies. However, we emphasize that some of the neurones that we recorded did not correlate with the activity of any of the muscles recorded. Such neurones might be related to controlling higher-level aspects of the movement or may even relate to more general kinematic variables (see Morrow et al., 2007).

The changes in phase and magnitude of cell discharge frequency of the overall population of task-related PTNs suffice to modify appropriately activity for those muscles in a synergy that are active in all three conditions (lead, control, and trail). A question arises, however, as to how muscles active only during the voluntary gait modifications are controlled. This is a particularly important consideration for those muscles that are strongly activated only in the trail condition, e.g., the Tri(1) burst. There are several possible ways in which this burst could be controlled. First, it is possible that cells that contribute to the production of activity in the extensor muscles during stance would discharge an additional burst of activity in the swing phase during the trail condition and contribute to activity both during swing and stance. However, an examination of the database of cells used in previous publications (Drew, 1993; Drew et al., 1996) showed little evidence for this pattern of activity. A second possibility is that there is a separate population of cells that discharge only, or primarily, in the trail condition and are responsible for the production of these additional bursts of activity. Again, there was little evidence for this in our population. A third possibility is that the discharge activity in the population of cells that influence activity in e.g., the TrM(1) is sub-threshold to produce activity in e.g., the Tri(1) during the unobstructed and lead conditions. The large increase in activity observed in these cells during trail, however, would be sufficient to produce a (non-linear) increase in the excitation of this muscle. Such a non-linear effect might be facilitated by a change in afferent input produced by the modified movement. It is interesting to note, for example, that low intensity stimulation of peripheral afferents during the swing phase of locomotion evokes strong responses in the Tri, despite the lack of any natural activity in the motoneurones at this time (Drew and Rossignol, 1987).

In addition to issues of the characteristics of the control signal, one also has to consider the anatomical bases for the proposed control system. If a given synergy is controlled as a unit by activity in a small homogeneous subpopulation of PTNs, then one would expect one or both of the following conditions to be met: (1) given that there are no monosynaptic connections from the motor cortex to motoneurones in the cat (Illert et al., 1976), PTNs should project to interneurons in the spinal cord that selectively activate the muscles identified as belonging to a single synergy (**Figure 10A**) and/or (2) cells in the motor cortex with more restricted projections to a subset of the muscles identified as belonging to a single synergy are connected by intra-cortical connections (**Figures 10B,C**). General evidence for both of these propositions exists. The anatomical studies of Alstermark (Tantisira et al., 1996) and the electrophysiological studies of Lundberg et al. (1987) emphasize that interneurons within the spinal cord branch to activate groups of proximal and distal muscles. In addition, Shinoda has very elegantly shown that individual corticospinal axons branch to multiple levels of the cervical spinal cord (Shinoda et al., 1976, 1986; Futami et al., 1979). Together, this provides two mechanisms by which small populations of cortical neurones could activate muscle synergies.

Interestingly, Hart and Giszter (2010) have used spike triggered averaging (STA) to demonstrate that the muscles activated by a given interneurone in the frog spinal cord correlate strongly with the weighting matrix of individual motor primitives, or synergies.

sequentially (filled histogram) as the cat steps over the obstacle. (Figure adapted from Drew et al., 2008b). **(B)** We illustrate activity of a PTN (different

> Furthermore, in the primate, STA demonstrates that individual corticomotoneuronal neurones in the motor cortex can influence the activity of multiple muscles including those acting at proximal and distal joints (McKiernan et al., 1998; Griffin et al., 2008).

magnitude and phase of the PTN and the TrM. The data for these illustrations

is taken from a previous publication (Drew, 1993).

Similarly, microstimulation of the motor cortex, particularly during locomotion (Armstrong and Drew, 1985) simultaneously activates both proximal and distal muscles. At the cortical level, there are abundant references to show the richness of the corticocortical connections between both adjacent and more distant regions of cortex (Huntley and Jones, 1991; Schneider et al., 2002; Capaday et al., 2009, 2011; Smith and Fetz, 2009). In addition, the motor cortex receives abundant input from premotor cortex and from the posterior parietal cortex (Ghosh, 1997; Andujar and Drew, 2007), which might also serve to coordinate activity between different subpopulations of PTNs.

Missing from the experimental data in the mammal, however, is the information as to whether the muscles linked anatomically or functionally in the studies mentioned above correspond to the synergies identified by our, or other, studies. In addition, there is no evidence of a direct link between the discharge pattern of individual cells and the muscles of the synergy that is activated or influenced by its discharge.

# **COMPARISON WITH OTHER MODELS OF MODULARITY**

Our model shares the same basic concept proposed by Bizzi (see Introduction for references), namely that a complex series of behaviors can be produced by the differential combination of a finite number of synergies. This base concept has been tested and confirmed in a large number of behaviors as listed in the Introduction. However, the majority of these models, whether based on synchronous or time-varying synergies (d'Avella et al., 2003), propose a small number of synergies, typically 4–6, although sometimes fewer, to control a wide range of complex movements. Moreover, the control signals, represented by the waveforms resulting from the mathematical decomposition, especially those obtained using principal component analysis (PCA), are generally quite broad and occupy a large proportion of the step cycle (e.g., Lacquaniti et al., 2012) or the movement phase (e.g., Overduin et al., 2008). These control signals activate muscles according to a weighting matrix. As such, all muscles are activated to varying amounts by all control signals and the movement produced is the sum result of the muscle activation although generally a small number of muscles has relatively large weights in each synergy (see e.g., Kargo and Nitz, 2003; Hart and Giszter, 2004). Moreover, the calculation of the synergies is generally performed by including as wide a range of movements as possible to ensure that the resultant synergies are capable of reproducing a large behavioral repertoire. It is generally considered that these signals are then transformed at the level of the spinal cord to produce the appropriate periods of muscle activity (see e.g., Ting, 2007 for a schema for postural control). Some support for this view of motor control comes from recent experiments by Overduin et al. (2012) who showed that the EMG responses evoked by long trains of stimulation (150–500 ms) at distributed points in the primate motor cortex could be reduced to three synergies. Moreover, the synergies evoked by microstimulation corresponded to the synergies defined by EMG activity during reach and grasp movements.

In contrast, our approach and resulting conceptual model differs in two important considerations from those of the more commonly used methods. The first is that we use a more classical method of defining a synergy (see above) and place no limit on the number of synergies that can be defined. During locomotion, this results in the identification of 10 or 11 synergies rather than the 4–6 more commonly used. Moreover, the large number of synergies and the fact that a single muscle burst cannot be included in more than one synergy results in the sparse synergies defined in this manuscript and in Krouchev et al. (2006). These synergies occupy only a small part of the step cycle and are compatible with a control mechanism by which each synergy is modulated by a specific sub population of cortical neurones as defined in the previous section. The second difference is that we define our synergies during only a single behavior, locomotion. This is based on the premise that the neural circuits in the spinal cord primarily evolved to be optimal for producing the basic pattern of locomotion. As such, we propose that the motor cortex (and other descending pathways) act via these circuits, adapting their activity to the behavior as required. In some cases, as demonstrated here for voluntary gait modifications, modifying the phase and magnitude of these synergies is sufficient to produce changes in behavior. In other behaviors, however, the challenge might well be to modify these synergies by changing both the number and the composition of these synergies (see also Kargo and Giszter, 2000, 2008). This contrasts with the idea of a limited number of immutable synergies that can be used to produce a wide range of behaviors (d'Avella and Bizzi, 2005; Overduin et al., 2008).

Our results can also be compared to other studies and in other tasks that have specifically studied the contribution of motor cortical discharge to the control of synergies. Kargo and Nitz (2003), for example, also showed motor cortical cells that were sequentially activated in rats trained to reach. They used independent component analysis (ICA) to identify a number of synergies active throughout a reach. Cross-correlation analyses showed that a majority of task-related motor cortical cells showed a significant correlation with only a single independent component (IC). When cross-correlated with the activity of individual muscles, neurones were significantly correlated with many, although not all of the muscles included in a synergy. Miller (Holdefer and Miller, 2002; Morrow et al., 2009) also suggested the existence of a number of muscle synergies during a reaching task, with each synergy controlled by separate populations of motor cortical neurones. These synergies were relatively stable across task (Morrow et al., 2009). In our own work (Yakovenko et al., 2011), using similar methods to those in this manuscript, we have also shown sparse synergies during reaching movements that are similar to those observed during locomotion and we have suggested that PTNs in the motor cortex control the limb trajectory during reaching in the same manner that we propose for voluntary gait modifications.

The model that we propose is clearly most similar to the original modular concepts proposed by Grillner (1981). As illustrated in Grillner and Wallen (1985), we propose that descending pathways act via modules in the spinal cord to modify locomotion. We further propose that the detailed motor commands observed in sub-populations of motor cortical neurones (Drew, 1993; Drew et al., 1996, 2004, 2008b) provide a means for exerting specific control over muscle activity while ensuring that the resulting gait modifications are fully integrated into the base locomotor rhythm (Drew, 1991a; **Figure 10D**). However, it should also be noted that in the majority of these experiments, again including our own, the synergies that are identified are multi-articular and not confined simply to a single joint as in the original model of Grillner (1981).

It could be argued that the synergies that we are defining in the intact cat emerge from the concerted activity of a unit CPG of the type produced by Griller, moulded by descending and peripheral afferents to produce our more fractionated, multi-joint modules. Indeed, a recent study by Markin et al. (2012), using methods similar to those described here, found some differences

**(A)** a single subpopulation of PTNs connects with spinal interneurons (ins) that project to the motoneurones (mns) of all muscles in a synergy. **(B)** Two subpopulations of PTNs, linked by intracortical connections, and discharging simultaneously during locomotion connect with different populations of spinal interneurons, each of which innervate only part of the total muscle synergy. **(C)** Same principal as **(B)** but providing the possibility of more fractionation. Dotted lines in **(A–C)** indicate weaker connections to interneurons and motoneurones that are not part of the synergy directly modified by the motor sub-populations of PTNs with spinal modules comprising the CPG. Each module will activate the motoneurones of a given synergy. Note that in this illustration we have separated the rhythm generator from the pattern generator elements as in the model by Rybak et al. (2006) and as supported by the results of our microstimulation studies (Rho et al., 1999). **(D)** is modified from Drew (1991a); Drew et al. (2008b). Abbreviations: E, extensor part of rhythm generator; F, flexor part of rhythm generator; in, interneurone; mn, motoneurone; PTN, pyramidal tract neurone; Sn, synergy #n.

in the synergies identified in the hindlimb of intact cats and those identified during fictive locomotion in decerebrate cats. These differences were primarily observed in muscles acting around two joints suggesting that the expression of the final pattern of activity in such muscles is modified by peripheral input. In this respect, it is important to realize that the patterns of activity observed in spinal animals are strongly modulated by peripheral input (Pearson and Rossignol, 1991; Rossignol et al., 1993; Lemay and Grill, 2004; Saltiel and Rossignol, 2004a,b; Cheung et al., 2005).

# **SYNERGIES AS A UNIT OF CONTROL OF GENERAL MOTOR ACTIVITY**

Even in cats a large number of forelimb movements do not involve the simple sagittal pattern of intralimb coordination that we have so far examined and in primates, and especially in humans, the patterns of limb movement become more and more flexible. This is especially true when we consider grasping movements that involve both arm movement and control of the hand. Indeed, the extent to which control of hand movements may be explained by synergies is quite controversial (Brochier et al., 2004; Theverapperuma et al., 2006; Kutch et al., 2008; Overduin et al., 2008). How then does the concept of a control of movement by synergies, at least those of the type that we propose, lend itself to a flexible control of limb movement throughout a wide behavioral range?

The answer to this question lies to some extent in whether one considers muscle synergies as a concept that simplifies motor control, as is generally assumed (see Tresch and Jarc, 2009) or one that is the result of an evolution in which more recently developed descending pathways must act through spinal circuits that have been conserved (see Krouchev et al., 2006; Giszter et al., 2007; Giszter and Hart, 2013). As mentioned above, we favor the second possibility. We suggest that with the parallel evolution of the nervous system and the musculoskeletal system, spinal circuits became adapted to provide more flexible and agile locomotor limb movements, culminating in reaching and ultimately, in primates, reaching and grasping movements (Georgopoulos and Grillner, 1989). In this view, as movements became more complex, the challenge for the nervous system was to produce flexible movements that are independent of the more stereotypical arm movements observed during locomotion. Indeed, from an evolutionary viewpoint one might speculate the development of a hierarchical control system that evolved together with the increasing flexibility required in the control of movement. Such a hierarchy is frequently discussed in the locomotor literature (Rossignol, 1996) but less frequently with respect to voluntary movements (see Ting, 2007; Roh et al., 2011).

During locomotion, it is well established that circuits in the spinal cord generate a rhythm that contains details about the pattern of muscle activity (see Introduction). These spinal circuits are then subject to modification by brainstem and cortical inputs. The brainstem pathways are suggested to regulate the level of muscle activity and posture, for example, during walking uphill (Orlovsky, 1972; Drew et al., 2004). This general control of muscle activity is facilitated by the fact that the axons of neurones in the brainstem pathways, the reticulo- and vestibulo-spinal tracts, have diffuse termination patterns that influence muscle activity around multiple joints, and in multiple

limbs (Matsuyama et al., 1988; Drew and Rossignol, 1990a,b). Moreover, the signals recorded from reticulospinal neurones during unobstructed locomotion (Drew et al., 1986; Matsuyama and Drew, 2000) and during voluntary gait modifications (Prentice and Drew, 2001) do not show the same level of fractionation as signals from the motor cortex. Indeed, many cells in the reticulospinal system show broad patterns of modulation that are more compatible with the waveforms identified by many of the decomposition studies than those observed in the motor cortex. Moreover, rather than specifying the exact patterns of motor activity that need to be produced, we have suggested that the signals from these pathways are integrated with the rhythmical signals in the spinal cord to modify the production of many of the more stereotypical patterns of behavioral activity, including locomotion (Drew et al., 2004). In some respects, therefore, the activity and connectivity of brainstem pathways resembles that postulated for synergies in general, and especially those implicating the brainstem (Roh et al., 2011). Nonetheless, it should be emphasized that there is little evidence for a limited subset of discharge patterns among reticulospinal cells during locomotion as might be expected if they represented a small range of synergies. Moreover, it should also be noted that synergy studies generally concentrate on the activity within a single limb whereas many neurones in the reticulospinal pathways appear to be involved in coordinating activity in two, or more, limbs.

At a higher level of the hierarchy, the motor cortex (and the red nucleus) provides a more specific level of control of locomotion. During unobstructed locomotion, the contribution is likely to be one of step-by-step regulation of the step cycle superimposed on the control exerted at lower levels of the nervous system. A contribution at this level is supported by the fact that motor cortical cells, including PTNs are modulated during unobstructed locomotion [reviewed in Drew et al. (1996)] and that microstimulation of the motor cortex can modify both the pattern and the timing of the step cycle (Armstrong and Drew, 1985; Rho et al., 1999). The major contribution of the motor cortex, however, is in the control of precise locomotor movements on the basis of visual information (Drew, 1988, 1993; Beloozerova and Sirota, 1993). As developed in this manuscript, we propose that the motor cortex produces gait modifications by altering the timing, duration and relative timing of muscle synergies. As suggested above, these signals are likely to be mediated by the same neuronal circuits responsible for the generation of the locomotor rhythm and the pattern of locomotion. This allows a specific control over the level of activity in the different synergies, thus adapting activity to the specific requirements of the task while still integrating this into the step cycle (**Figure 10D**). In our model, control of voluntary gait modifications does not require the production of new synergies but rather modification of the synergies already present in unobstructed locomotion. This should be compared with the results of Ivanenko et al. (2005) in which stepping over an obstacle by human subjects required the addition of a new principal component that was responsible for explaining >20% of the overall variance in the EMG patterns.

At the highest level there is a need to control movements that are non-stereotypic. This might involve breaking apart synergies or the addition of new synergies. For example, during corrective manoeuvres, Kargo and Giszter (2000) have shown there is a need to modify synergy composition to correct for perturbations. In other conditions, for example when making relatively simple movements, there might be a need to activate 1 or 2 sparse synergies independently of others. Other movements may require separation of synergies, for example, uncoupling activity in wrist and elbow muscles. This may be achieved by activating only a part of the population of PTNs contributing to muscle activity in a given synergy (e.g., **Figure 10B**). Additionally, it might require inhibiting some of the motoneurones that would normally be activated by a given sub-population of PTNs. Indeed, it has been suggested that the development of the motor cortex and, in primates, the corticomotoneuronal system, provides the ability to generate more fractionated movements, in part by inhibiting some synergies and activating others (Lemon and Griffiths, 2005; Drew et al., 2008a). In other movements, particularly those involving the fingers, the number of patterns of motor cortical recruitment makes it difficult to argue that movement is controlled by synergies. Poliakov and Schieber (1999), for example, failed to find any evidence of broad groups of motor cortical neurones controlling finger movements and instead reported a highly diverse pattern of activity suggesting control by a distributed network.

# **REFERENCES**

*J. Neurophysiol.* 92, 1770–1782. doi: 10.1152/jn.00976.2003

# **SUMMARY**

The data illustrated in this manuscript are compatible with a view that motor cortical cells contribute to the control of locomotion by modulating the activity of a limited number of sparse synergies. We suggest that a relatively simple control system can serve to control limb trajectory under a wide range of situations simply by modifying the phase and magnitude of the period of activity in a limited number of functionally segregated populations of PTNs. In addition, Yakovenko et al. (2011) showed that similar synergies are observed during reach as during locomotion. Indeed, the sequential organization of sparse synergies in unobstructed locomotion, stepping over obstacles and even reaching, is strongly suggestive of a relationship to prime mover muscles and their agonists. These relationships, in turn, are activated subject to biomechanical (kinematic, dynamic) constraints and to existing neural circuits provided by subcortical and spinal neural patternand rhythm- generating circuitry.

# **ACKNOWLEDGMENTS**

This work was supported by a CIHR operating grant (MOP 9558) to T. Drew, a New Emerging Team Grant in Computational Neuroscience from the CIHR, and an infrastructure grant from the FRSQ. We thank Drs. Elaine Chapman and John Kalaska for their comments on this manuscript.

muscle synergies during natural muscle behaviors. *J. Neurosci.* 25, 6419–6434. doi: 10.1523/ JNEUROSCI.4904-04.2005

synergies in the construction of a natural motor behavior. *Nat. Neurosci.* 6, 300–308. doi: 10.1038/ nn1010

165, 323–346. doi: 10.1016/S0079- 6123(06)65020-6

muscle synergies during locomotion in the intact cat as revealed by cluster analysis and direct decomposition. *J. Neurophysiol.* 96, 1991–2010. doi: 10.1152/jn.00241.2006

linkage of reflex actions to alphamotoneurones. *Exp. Brain Res.* 65, 271–281. doi: 10.1007/BF00236299

*Neuron* 76, 1071–1077. doi: 10.1016/j.neuron.2012.10.018

circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. *J. Physiol.* 577, 617–639. doi: 10.1113/jphysiol.2006.118703

the isolated mouse spinal cord during spontaneous deletions in fictive locomotion: insights into locomotor central pattern generation organization. *J. Physiol.* 590, 4737–4759. doi: 10.1113/jphysiol.2012.240895

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 26 March 2013; paper pending published: 16 April 2013; accepted: 12 June 2013; published online: 11 July 2013.*

*Citation: Krouchev N and Drew T (2013) Motor cortical regulation of sparse synergies provides a framework for the flexible control of precision walking.* *Front. Comput. Neurosci. 7:83. doi: 10.3389/fncom.2013.00083 Copyright © 2013 Krouchev and Drew. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.*

# Increased motor cortex excitability during motor imagery in brain-computer interface trained subjects

### *Olesya A. Mokienko1,2, Alexander V. Chervyakov1 \*, Sofia N. Kulikova1, Pavel D. Bobrov2,3, Liudmila A. Chernikova1, Alexander A. Frolov2,3 and Mikhail A. Piradov1*

*1 Research Center of Neurology Russian Academy of Medical Science, Moscow, Russia*

*2 Institute of Higher Nervous Activity and Neurophysiology of RAS, Moscow, Russia*

*3 Technical University of Ostrava, Ostrava, Czech Republic*

### *Edited by:*

*Yuri P. Ivanenko, IRCCS Fondazione Santa Lucia, Italy*

### *Reviewed by:*

*Mikhail Lebedev, Duke University, USA*

*Francesca Sylos Labini, IRCCS Santa Lucia Foundation, Italy*

### *\*Correspondence:*

*Alexander V. Chervyakov, Research Center of Neurology Russian Academy of Medical Science, Marshala Katukova Street, 13-2-124, 123181 Moscow, Russia e-mail<EMAIL_ADDRESS>
**Background**: Motor imagery (MI) is the mental performance of movement without muscle activity. It is generally accepted that MI and motor performance have similar physiological mechanisms.

**Purpose**: To investigate the activity and excitability of cortical motor areas during MI in subjects who were previously trained with an MI-based brain-computer interface (BCI).

**Subjects and Methods**: Eleven healthy volunteers without neurological impairments (mean age, 36 years; range: 24–68 years) were either trained with an MI-based BCI (BCI-trained, *n* = 5) or received no BCI training (*n* = 6, controls). Subjects imagined grasping in a blocked paradigm task with alternating rest and task periods. For evaluating the activity and excitability of cortical motor areas we used functional MRI and navigated transcranial magnetic stimulation (nTMS).

**Results:** fMRI revealed activation in Brodmann areas 3 and 6, the cerebellum, and the thalamus during MI in all subjects. The primary motor cortex was activated only in BCI-trained subjects. The associative zones of activation were larger in non-trained subjects. During MI, motor evoked potentials recorded from two of the three targeted muscles were significantly higher only in BCI-trained subjects. The motor threshold decreased (median = 17%) during MI, which was also observed only in BCI-trained subjects.

**Conclusion**: Previous BCI training increased motor cortex excitability during MI. These data may help to improve BCI applications, including rehabilitation of patients with cerebral palsy.

**Keywords: brain-computer interface, motor imagery, navigated TMS, functional MRI, neurorehabilitation**

# **HIGHLIGHTS**

# **INTRODUCTION**

Modularity is important concept in understanding mechanisms of motor control and motor learning. Recent investigations of muscle synergies, motor primitives, compositionality, basic action concepts, and related work in machine learning have contributed to advance our understanding of the architecture underlying rich motor behaviors (Lacquaniti et al., 2013). One of the most interesting topic in researches of motor system modularity and organization is study of motor imagery and changes in neural networks during it.

 These structures include the premotor and supplementary motor cortices (Brodmann area 6), parietal cortical areas, cingulate gyrus, basal ganglia, and cerebellum. The primary motor cortex (Brodmann area 4) is also active during MI (Jeannerod, 2001; Neuper et al., 2005; Sharma et al., 2006; Dickstein and Deutsch, 2007; Mulder, 2007). Furthermore, several studies using transcranial magnetic stimulation (TMS) demonstrated increased corticospinal excitability and increased amplitudes of motor evoked potentials (MEPs) during MI (Fadiga et al., 1999; Hashimoto and Rothwell, 1999; Vargas et al., 2004; Cicinelli et al., 2006; Stinear et al., 2006; Pichiorri et al., 2011).

These previous findings led scientists to develop an MI training paradigm to stimulate neuroplastic changes in patients with paresis resulting from brain injury, or for use in athletic training. An electroencephalography (EEG)-based brain-computer interface (BCI) is a promising method to support MI during such training. The BCI transforms EEG signals generated during MI into commands that can control an external device (Prasad et al., 2010; Mokienko and Chernikova, 2011; Shih et al., 2012). Modulation of the sensorimotor rhythm can serve as the signal of brain activity during MI (Pfurtscheller and Lopes da Silva, 1999). However, the plasticity-related changes resulting from BCI-supported MI training have not been studied in detail. Furthermore, no study has included the navigated transcranial magnetic stimulation (nTMS) method, or compared functional magnetic resonance imaging (fMRI) and nTMS data for functional mapping during MI in BCI-trained and not trained subjects. The aim of our experiment was to investigate the activity and excitability of different cortical motor areas during MI in BCI-trained and BCI-naïve subjects.

# **MATERIALS AND METHODS PARTICIPANTS**

The inclusion criteria were an age of 20–70 years old, an absence of neurological disorders, normal vision, right-handedness, and written informed consent. Eleven volunteers (mean age, 36 years; age range, 24–68 years; 7 males and 4 females) were included into the study. Subjects of group 1 (*n* = 5, mean age = 45.8 y.o) had 10 to 15 sessions of BCI-supported training 20–30 min each. The training course was followed by fMRI and nTMS examinations. Subjects of group 2 (*n* = 6, mean age = 27.6 y.o) were tested without preliminary training session. All subjects underwent fMRI and nTMS after the training sessions. The protocol was approved by the local Ethical Committee of the Research Center of Neurology of RAMS, Moscow. All subjects provided written informed consent.

# **BCI TRAINING**

The BCI training was based on EEG activity patterns recording during grasping MI. Subjects sat comfortably in an armchair located 1 m from a computer screen that presented visual instructions. Subjects visually fixated on a circle presented in the center of the screen and received instructions from three surrounding figures (rhomboidal arrows). Subjects were given three commands instructing them to relax (upper arrows were illuminated) or imagine slow grasping movements with the right or left hand (right or left arrow illuminated). The "Relax" command meant that the subject had to sit still and look at the center of the screen. Commands were presented randomly, each of 10-s duration. For each subject, training was performed in 10–15 experimental days, with one 20–30 min session performed each day. Intervals between training sessions were 1–4 days.

A visual cue provided the subject with feedback regarding the mental task recognition: the central circle turned green if the classifier recognized the task in agreement with the given command, or remained white if the signal was not recognized. The EEG was registered with 30 electrodes distributed over the head in accordance with the standard international 10–20 system. EEG signals were filtered from 5–30 Hz. We used a Bayesian approach for EEG pattern classifying. The activity sources most relevant for BCI functioning were identified using an independent component analysis (ICA). Classification accuracy was measured with Cohen's kappa, a parameter conventionally used in BCI studies (Kohavi and Provost, 1998). A kappa of 1 indicates perfect recognition, whereas a kappa of 0 indicates random recognition.

# **fMRI**

fMRI was conducted using a Magnetom Avanto 1.5-T MRI system (Siemens, Germany). Standard axial T2-weighted turbo-spin echo imaging was performed initially to rule out pathological changes in brain tissue [repetition time (TR), 4000 ms; echo time (TE), 106 ms; section thickness, 5.0 mm; matrix, 230 × 230 mm; imaging time, 2 min 2 s]. Anatomical data were obtained with sagittal T1-weighted gradient echo imaging with isometric voxels (Ò1 multiplanar reconstruction: TR, 1940 ms; TE, 3.1 ms; TI, 1100 ms; section thickness, 1.0 mm; matrix, 256 × 256 mm; imaging time, 4 min 23 s). During the fMRI experiment, subjects performed the same task that was performed during the BCI training sessions, but without feedback. For each subject, three sets of functional data were obtained representing different conditions, including rest (8 replicates) and right or left hand movement imagery (4 replicates each). The imaging mode used was axial T2∗ gradient echo (TR, 3800 ms; TE, 50 ms; matrix, 192 × 192 mm, section thickness, 3 mm) with fat suppression and correction for motion. The imaging time was 6 min 10 s.

The data analysis was performed in the MATLAB (Mathworks, Natick, MA, USA) environment using the statistical package for processing in SPM8 (Welcome Trust Center of Neuroimaging, London, UK). The first step of analysis corrected head movement artifacts. Next, the functional data were translated to Montreal Neurological Institute (MNI) coordinates (i.e., normalization). Standard MNI coordinates, which are used in the SPM8 package, were developed and adopted by the International Consortium for Brain Mapping. In the next step, the normalized data were smoothed. This step was followed by a classic analysis that used generalized linear models. The results from each subject were used in the group analysis to identify areas showing task-specific activity.

# **nTMS**

Neurophysiological investigation was performed using nTMS with the NBS eXimia Nexstim apparatus (Finland).It includes 70 mm figure-eight-shaped BiPulse Nexstim coil, with a maximal magnetic field strength of 199 V/m and a magnetic impulse duration of 280µs. The coil was placed anteromedially at a 45◦ angle from the midline. The stimulated hemisphere was not the same for all the subjects and was chosen randomly.

As a first step, all subjects underwent an MRI investigation on a Magnetom Symphony 1.5 T scanner (Siemens, Germany) using a T1 multiplanar reconstruction regime (MPR); the data were loaded into the NBS eXimia Nexstim system to obtain subjects' individual 3D brain models. Following that, real anatomical entities were matched to their MRI representations.

The MEPs were recorded using a standard EMG machine (Nexstim, EMD, Finland) and surface electrodes. MEPs were recorded by placing 0.6 cm2 EMG electrodes on the target muscle being mapped [abductor pollicis brevis (APB), flexor carpi ulnaris (FCU), and extensor carpi radialis (ECR)] which were positioned according to the atlas of Leis and Trapani (2000), according to the belly–tendon principle. The ground electrode was placed on the right clavicle or on the upper third of the right forearm. We then determined the resting motor threshold (MT), defined as the lowest stimulation intensity allowing evocation of motor responses 0.50 mV peak to peak amplitude in 5/10 trials with the patient at rest (Rossini et al., 1994). Resting MT was measured in present (%) of the maximum intensity of the magnetic stimulator (1,5 Tesla). Evoked motor responses (EMRs) and their amplitudes and latencies were recorded for each target muscle. Cortical motor representations were constructed from these observations.

In the first step of the experiment, the areas of interest (contralateral primary motor and premotor cortices) were stimulated with magnetic fields of 80–110 V/m to identify EMRs with amplitudes of 100–500µV. The resting MT was determined for each site by the largest detected EMR amplitude. Cortical representations of the target muscles were mapped at 110% intensity of the determined resting MT. The mean EMR amplitudes and muscle motor representations were evaluated during cortical mapping.

In the second step of the experiment, the passive EMR threshold was determined and the motor representations were mapped while subjects imagined grasping with the contralateral hand. The motor representations were mapped using stimulus intensities as in the first step. The hand was positioned on the armrest in the neutral position of the radiocarpal joint.

Statistical analysis of quantifiable data was performed using a repeated measures analysis of variance (ANOVA) and Newman-Keuls *post-hoc* test using the Statistica 6.0 software package (StatSoft, 2003). The data are presented as the median and 25–75% quartiles. Differences were considered significant at *p <* 0*.*05.

# **RESULTS**

### **SUBJECTS**

Eleven volunteers (mean age, 36 years; age range, 24–68 years; 7 men, 4 women) participated in the study. Subjects in group 1 (*n* = 5) underwent 10 to 15 sessions of BCI-supported training that were each 20–30 min in duration. The training course was followed by fMRI and nTMS examinations. Subjects in group 2 (*n* = 6) were tested without performing preliminary training sessions.

### **BCI TRAINING**

The achieved accuracy rates (median Cohen's kappa) were 0.46 [0.45; 0.52]. BCI control for all subjects was achieved with sensorimotor rhythm modulation. MI was accompanied with desynchronization of mu and low beta rhythms (i.e., eventrelated desynchronization) (Pfurtscheller and Lopes da Silva, 1999). This signal was a recognizable command for the BCI.

### **FUNCTIONAL MRI**

In BCI-trained subjects, MI was accompanied by activity in the contralateral somatosensory (Brodmann area 3), primary motor (Brodmann area 4), and premotor cortical areas. Activity also was observed in the bilateral supplementary motor cortex (Brodmann area 6), contralateral ventral lateral nucleus of thalamus, and ipsilateral cerebellum (*p <* 0*.*0005; **Figure 1**).

For BCI-naïve subjects, activity was observed in the contralateral somatosensory (Brodmann area 3) and premotor cortical areas, as well as the supplementary motor cortex (Brodmann area 6) bilaterally. Other activated areas included the contralateral ventral lateral nucleus of the thalamus, ipsilateral cerebellum, contralateral Brodmann area 9, and bilateral Brodmann areas 40 and 13 (*p <* 0*.*0005). The primary motor cortex was not activated in untrained subjects (**Figure 2**). The areas of activation including the somatosensory, premotor, and supplementary motor areas observed during MI were significantly larger in BCI-naïve subjects than in BCI-trained subjects (*p <* 0*.*01).

### **nTMS**

In BCI-trained subjects, the passive MT for the motor cortex decreased by 6–18% (median change was 17%) during MI compared to the rest condition. In BCI-naïve subjects, the threshold change during MI compared to rest were not significant and were inconsistent among subjects. The threshold decreased by 1–8% in three subjects, insignificantly increased in two subjects, and was unchanged in one subject. The MT changes were statistically significant only for BCI-trained subjects (*p <* 0*.*01, **Table 1**).

For APB, the median change in motor response during MI compared to rest condition was 63% in BCI-trained subjects, and 11% in BCI-naïve subjects. For ECR, the change was 150% in BCI-trained subjects and 1% in BCI-naïve subjects. In BCItrained subjects, the responses in APB and ECR (mean EMR) during MI were significantly higher during MI compared to the rest condition (APB, *p* = 0*.*03; ECR, *p* = 0*.*01). In BCI-naïve subjects, the differences in EMR were not significant (APB, *p* = 0*.*24; ECR, *p* = 0*.*23, **Table 1**).

We did not observe statistically significant increases in mean EMR amplitude in FCU for either group (**Table 1**). The median change in motor response was 78% in BCI-trained subjects

**FIGURE 2 | Areas of activation during grasping imagery in untrained subjects (group analysis of fMRI data, Left hand imagery** *>* **Rest,** *p <* **0***.***0005). (A),** Brodmann areas 3 and 4; **(B)**, supplementary motor cortex; **(C)**, cerebellum; **(D)**, insula; **(E)**, Brodmann area 9; **(F)**, Brodmann area 40.

**Table 1 | Motor thresholds and evoked motor responses during rest and motor imagery for the two groups (represented as median, [25th and 75th percentiles]).**

*EMR, evoked motor response; APB, abductor pollicis brevis; FCU, flexor carpi ulnaris; ECR, extensor carpi radialis. ANOVA with Newman-Keuls post-hoc test.*

and 12% in BCI-naïve subjects. Moreover, in BCI-trained subjects, stimulation induced EMRs that were larger during MI than at rest, which was not observed in BCI-naïve subjects (**Figure 3**). A comparison of nTMS and fMRI maps revealed partial overlap of motor areas detected by these two methods (**Figure 4**).

# **DISCUSSION**

The activation and excitability of the motor cortex during MI is different in BCI-trained and BCI-naïve subjects, and this difference can be detected with a combination of fMRI and nTMS.

# **MI AND ACTIVATION OF MOTOR STRUCTURES**

The fMRI analysis revealed the same brain areas were active during MI for both groups. These areas included the contralateral somatosensory, contralateral premotor, supplementary motor cortex bilaterally, contralateral ventral lateral nucleus of thalamus, and ipsilateral cerebellum. Similar activation patterns were described in other fMRI-based MI studies (Jeannerod, 2001; Neuper et al., 2005; Sharma et al., 2006; Dickstein and Deutsch, 2007; Mulder, 2007).

In the literature, there is a debate regarding what role the primary motor cortex plays in MI, as some studies failed to observe its activation (Parsons et al., 1995; Hanakawa et al., 2003; Meister et al., 2004; de Lange et al., 2005). In our study, primary motor cortex activation was observed only in BCI-trained subjects. Therefore, we suppose primary motor cortex is involved in MI in individuals who can "successfully" imagine a movement, or who have been trained to do so (e.g., BCI-supported training).

The somatosensory, premotor, and supplementary motor cortical areas were activated during MI, and were larger on fMRI in BCI-naïve subjects. This is in agreement with the principles implying localization of new or skilled movements.

In BCI-naïve subjects, we observed bilateral activation of Brodmann area 40, the complex associative cortex. This region plays a central role in developing cognitive strategies and motor programs, and its activation was described in several previous MI studies (Gerardin et al., 2000; Lafleur et al., 2002; Jackson

et al., 2003). This associative area was reported to be activated predominantly in the left hemisphere during complex motor performed by right-handed individuals (Gerardin et al., 2000). In addition, the dorsolateral prefrontal cortex (Brodmann area 9) was active in BCI-naïve subjects. This associative area represents the highest level of motor planning and regulation, and plays an important role in sensory and mnemonic information integration and working memory processes (Derrfuss et al., 2004). Right and left insula activation can be associated with cognitive control, task coordination, and working memory involvement (Derrfuss et al., 2004). It should be noted that BCI-trained subjects did not show significant activity in the associative areas.

# **MI AND CORTICOSPINAL EXCITABILITY**

Our nTMS findings indicate that MI is generally associated with a decrease in the evoked response threshold, an increase in the EMR amplitude, and an expansion of evoked response areas against the background of decreased excitation thresholds. Together, these changes reflect increased motor cortex excitability during MI. These changes are increased with MI training and often do not occur in untrained individuals. Our results are in agreement with the findings of other studies using classical TMS (without MRI navigation) (Fadiga et al., 1999; Hashimoto and Rothwell, 1999; Vargas et al., 2004; Cicinelli et al., 2006; Stinear et al., 2006).

### **nTMS AND INVESTIGATION OF MI**

Pichiorri et al. (2011) used TMS to assess the neuroplastic changes associated with MI-based BCI training. In that study, 10 healthy volunteers participated in 6–8 40-min BCI sessions. The training resulted in a significant increase in motor cortex excitability, and enhanced EMRs in target muscles during MI (Pichiorri et al., 2011). The TMS used in those studies was not navigated using MRI or fMRI data. In contrast to conventional TMS, nTMS allows local and precise stimulation based on an individual's MRI data (Chervyakov et al., 2013). This technique makes it possible to assess cortical excitability with a high spatial (2 mm) and temporal resolution. In the present study, we obtained similar results in terms of EMR, but we used both nTMS and fMRI. nTMS allowed us to map target muscle representations during MI for each subject based on MRI and fMRI data.

nTMS can be used to evaluate the dynamics of neuroplastic processes accompanying MI. MI mapping is most commonly performed with fMRI. In this case, the indirect measure of brain

**FIGURE 4 | Comparison of mapping results obtained during motor imagery in a BCI-trained subject by using fMRI (fist clenching imagery task) and nTMS.**