Correlation between cortical activity and electromyographic (EMG) activity of limb muscles has long been a subject of neurophysiological studies, especially in terms of corticospinal connectivity. Interest in this issue has recently increased due to the development of brain-machine interfaces with output signals that mimic muscle force. For this study, three monkeys were implanted with multielectrode arrays in multiple cortical areas. One monkey performed self-timed touch pad presses, whereas the other two executed arm reaching movements. We analyzed the dynamic relationship between cortical neuronal activity and arm EMGs using a joint cross-correlation (JCC) analysis that evaluated trial-by-trial correlation as a function of time intervals within a trial. JCCs revealed transient correlations between the EMGs of multiple muscles and neural activity in motor, premotor and somatosensory cortical areas. Matching results were obtained using spike-triggered averages corrected by subtracting trial-shuffled data. Compared with spike-triggered averages, JCCs more readily revealed dynamic changes in cortico-EMG correlations. JCCs showed that correlation peaks often sharpened around movement times and broadened during delay intervals. Furthermore, JCC patterns were directionally selective for the arm-reaching task. We propose that such highly dynamic, task-dependent and distributed relationships between cortical activity and EMGs should be taken into consideration for future brain-machine interfaces that generate EMG-like signals.
- brain-machine interface
- motor cortex
- neural dynamics
despite many decades of research, the way that cortical signals control voluntary muscle contractions is not well understood (Kalaska 2009; Levine et al. 2012; Rothwell 2012; Schieber 2011). Particularly, it is still controversial (Hirasaki et al. 2004) whether motor cortical neurons encode muscle force (Cheney and Fetz 1980; Kalaska 2009; Levine et al. 2012; Rothwell 2012), limb kinematics (Georgopoulos et al. 1982; Moran and Schwartz 1999; Oby et al. 2013) or higher-order representations of movements (Ethier et al. 2012; Wessberg et al. 2000;).
Spike-triggered averaging (STA) of EMGs has been an influential approach for several decades in the analysis of functional connectivity between cortical neurons and spinal motoneurons (Fetz et al. 1976; Fetz and Cheney 1980; Mantel and Lemon 1987; Muir and Lemon 1983). This approach has proven to be useful to detect short-latency mono- and polysynaptic connections between cortical neurons and spinal motoneurons during different motor behaviors (Baker et al. 1997; Buys et al. 1986; Davidson et al. 2007; Jackson et al. 2003; McKiernan et al. 2000; Muir and Lemon 1983; Smith and Fetz 1989). STA usually represents average comodulations taken from many different behavioral epochs (McKiernan et al. 1998; Semmler et al. 2000), even though STA variations between different epochs of behavioral tasks (Baker et al. 1997; Morrow and Miller 2003) and different tasks (Zhou and Fuster 2000) have also been reported.
Typically, STA is used to detect influences of cortical motoneurons on their target muscles during voluntary movements (Buys et al. 1986; Lemon et al. 1986), and to detect such influences before overt movements begin (Fetz et al. 2002). The same approach has also shown prominent correlations between somatosensory cortical neuronal firing modulations and muscle activity (Widener and Cheney 1997). Additionally, STA and cross-correlation analysis have revealed task-dependent entrainment of EMG activity to oscillations of cortical field potentials in rodents (Carvell et al. 1996; Nicolelis et al. 1995), nonhuman primates (Baker et al. 1997) and humans (Hari and Salenius 1999; Kilner et al. 2000; Mima et al. 2001; Mima and Hallett 1999).
Recently, the relationship between the activity of motor cortical neurons and EMGs came into the scope of research on brain-machine interfaces (BMIs). BMIs strive to reconstruct limb movements from the collective extracellular activity of simultaneously recorded cortical neuronal populations (Chapin et al. 1999; Nicolelis 2001; Wessberg et al. 2000). Several BMI studies have also focused on the reconstruction of muscle activity from neural ensemble recordings (Fitzsimmons et al. 2009; Oby et al. 2013; Santucci et al. 2005; Stark and Abeles 2007; Westwick et al. 2006). Overall, a variety of decoding algorithms have been used to extract EMGs from cortical signals. Linear decoders are common, due to their simplicity and ease of use, as well as their effectiveness in extracting behavioral variables from neuronal activity (Lebedev et al. 2008; Wessberg et al. 2000). Linear algorithms applied to decoding EMG signals from cortical activity include the Wiener filter (Carmena et al. 2003; Fitzsimmons et al. 2009; Santucci et al. 2005) and multiple input-single output systems (Pohlmeyer et al. 2007), both of which make predictions based on a combination of weighted neural activity (Pohlmeyer et al. 2007). Other approaches incorporate nonlinear elements, such as the Wiener cascade, which is a dynamic linear system that includes a fixed nonlinear element (Pohlmeyer et al. 2009). While the goal of this research is to produce the best transfer function from neuronal activity to EMGs, the exact relationship between individual neurons and EMGs is often buried in the decoder parameters and is difficult to visualize or analyze.
The present study adds several novel findings to the existing literature on cortico-EMG correlations and BMIs that mimic muscle force. We started with a hypothesis that correlations between the activity of cortical neurons and EMGs of arm muscles could change on a timescale of 50–100 ms, i.e., a typical interval for changes to occur in motor tasks executed both naturally and through a BMI. With the exception of one study on sinusoidal tracking movements using cross-correlation techniques (Schwartz and Adams 1995), previous studies of cortico-EMG relationships did not consider such short timescales. These studies also rarely related correlation patterns to specific task events, such as instructed delays and movement onset, or to spatial parameters, such as movement direction. Here, we utilized a technique we call joint cross-correlation (JCC) analysis to address these questions. This method was adapted from the joint peristimulus time histogram (JPSTH) analysis, a method originally introduced to measure correlations between pairs of neurons (Brown et al. 2004; Cardoso de Oliveira et al. 2001; Gerstein and Perkel 1972; Palm et al. 1988), but never applied to cortico-EMG correlations.
In the present study, multisite neural ensemble recordings were obtained using multielectrode arrays chronically implanted in multiple cortical areas, including the primary motor cortex (M1), the dorsal premotor cortex (PMd), and the primary somatosensory cortex (S1). The behavioral task consisted of either self-timed touch pad presses or center-out reaches with a joystick. JCC was employed to detect dynamic changes in the correlation between surface recordings of EMGs of multiple muscles and the activity of ensembles of cortical neurons. JCCs determined these changes with a 50-ms time resolution, which is cruder than STA resolution for short-latency postspike effects in intramuscular EMGs (Fetz et al. 1976; Fetz and Cheney 1980; Mantel and Lemon 1987; Muir and Lemon 1983), but still sufficient for the detection of time-dependent changes in cortico-EMG correlation.
Using this approach, we found that cortical ensembles were correlated with EMGs of multiple muscles. The patterns of correlations between the neurons and EMGs were strongly time dependent and also directionally dependent in the reaching task. Matching results were obtained with STA methods applied to the same data. However, the STA approach was less useful for detecting time dependency of correlations because its time resolution for detecting changes was lowered by a relatively wide 500-ms data window. Additionally, as revealed by an analysis of shuffled behavioral trials, STAs reflected common stimulus-related modulations, rather than purely effective connectivity (Aertsen et al. 1989). To correct for this effect, we subtracted STAs calculated with shuffled trial data from the original STAs. Still, even these shuffle-corrected STAs were less accurate than JCCs in terms of depicting time dependency of neuronal-EMG correlations.
Overall, our results support a highly distributed and highly dynamic correspondence between cortical neuronal ensembles and spinal motoneurons. While the underlying neuronal circuitry is difficult to determine from these data alone, our analytic approach provides a practical tool for examination of input-output relationships for BMIs that intend to mimic EMG activity.
All experimental procedures were approved by the Institutional Animal Care and Use Committee (IACUC) and were conducted in accordance with the National Research Council's Guide for the Care and Use of Laboratory Animals (2011).
Implants and electrophysiological recordings.
The experiments were conducted in three rhesus macaques (Macaca mulatta; one female, two males) implanted with cortical multielectrode arrays (Fig. 1). Following the IACUC-approved protocol, animals were pair-housed before implantation and then later housed in separated home cages (44 × 33.5 × 33.5 in. in length × width × height, with an additional 2.5 in. in height for a waste pan). Monkey A was implanted in the primary motor area (M1) and PMd of both hemispheres (Fig. 1B). Each of the four electrode arrays was composed of 23 tungsten microwire electrodes, 35 μm in diameter. The microelectrodes were in either a rectangular or circular arrangement with 1 mm of separation between adjacent electrodes. Monkey B was implanted in M1, PMd, supplementary motor, and posterior parietal areas of the right hemisphere, while monkey C was implanted in M1 and primary somatosensory area (S1) in both hemispheres. Monkey B was implanted with four 32-electrode stainless steel microwire arrays grouped into 16 pairs. The electrodes were 45 μm in diameter. We only used recordings from the M1 array in monkey B for the present study (Fig. 1C). Monkey C was implanted with four 96-electrode stainless steel microwire arrays: one in the arm representation, and one in the leg representation of both hemispheres. Electrodes were arranged in two 4 × 4 grids of electrode triplets staggered at 300-μm intervals with one grid over M1 and the other over S1. Each electrode was 50 or 63.5 μm in diameter. We only used recordings from neurons in the right hemisphere in monkey C for this study (Fig. 1D).
To access the cortical areas, craniotomies were performed over the cortical areas of interest under inhalational anesthesia induced by isoflurane (Oliveira and Dimitrov 2008). The location of the craniotomy was determined using stereotaxic guidance and cortical surface landmarks, which included the central sulcus, arcuate sulcus and superior precentral dimple. Bone screws were inserted into the skull to ensure fixation of the implants. After implantation, microelectrode arrays were glued to the skull with dental cement. The microwire arrays, their connectors (Omnetics, Minneapolis, MN) and the implant assembly were affixed with dental cement. In monkey C, the implant was protected by a custom plastic headcap in the shape of an inverted cone covered by a lid.
Concurrent neuronal ensemble recordings were conducted using a 128-channel multichannel acquisition processor (MAP; Plexon, TX). Neuronal units were sorted using the MAP software according to shape templates and voltage thresholding. To assess whether recorded units were single or multiunits, we estimated refractory period using interspike intervals. Units that exhibited refractory periods greater than 1.6 ms (Fitzsimmons et al. 2009; Hatsopoulos et al. 2004) were considered single units. According to this criteria, 78.3% of units in monkey A were single units, and 21.7% were multiunits. In monkey B, 78.9% were single units, and 21.1% were multiunits. In monkey C, 57.5% were single units, and 42.5% were multiunits. All example units in Figs. 2, 3B, 5, 6, 9, A and B, and 10 were categorized in this way as single units and will be referred to as such or as neurons. Population plots included data from both single and multiunits, so recorded units in these plots will be referred to simply as “neural units” or “units.”
EMG signals were recorded with surface electrodes (DE-2.1, Delsys, Boston, MA) from the biceps (long head), triceps (long), wrist flexors (flexor carpi radialis and adjacent muscles) and wrist extensors (extensor carpi radialis longus and adjacent muscles) of the working arm in all monkeys, but data from the biceps in monkey C were not used due to poor signal quality. EMGs were amplified (10,000 × gain) and band-pass filtered (20–450 Hz) with a Delsys amplifier.
Each monkey was seated in a primate chair during behavioral tasks. Light and flexible cables connected their implants to head stage amplifiers and then to the MAP. The monkey's head was not restrained for recordings in monkeys A and C and was restrained by attaching the implant to an experimental frame in monkey B.
Monkey A was trained to perform the self-timed touch pad presses (Fig. 2). The data collected from this monkey was used in one previous publication (Lebedev et al. 2008), with no overlap in the reported results with the present study. During training and recording sessions, the monkey's left arm was restrained with Velcro bands, and the right shoulder was loosely restrained to allow for performance of the task with the hand. The monkey was trained to press and hold its hand over a touch pad (an aluminum cylinder that served as a touch sensor based on electrical resistance), then release after an allotted time. The force with which the monkey pressed on the touch pad was not measured. The touch pad was located at the monkey's waist level. To obtain a juice reward, the monkey was required to hold its right hand over the touch pad for anywhere between 2.5 s and 4.5 s and then release. If this task was accomplished, the monkey received a reward at the time of release. The monkey received no cues for when it should release other than the juice reward for correct performance.
Monkeys B and C were trained to perform center-out reaching movements with a two-degree-of-freedom joystick (Fig. 3). The joystick was positioned at the monkeys' waist levels, and both monkeys held it with their left hands. The monkeys were required to exert 8–9 N of force on the end of the joystick to move it through the maximum range of motion (30° from vertical). An infrared sensor on the top of the joystick indicated whether or not the monkey's hand was on the joystick. A 47.5 cm × 27 cm LCD monitor was located 55–65 cm in front of the monkey. It displayed a circular cursor (ring of diameter 0.5 cm) and a larger circle (7 cm in diameter) that served as targets of reaching movements. Forward and backward movements of the joystick relative to the monkey translated to up and down movements of the cursor, respectively. The monkeys started each trial by holding the cursor in a central target for 500–2,000 ms. The central target was then replaced by a peripheral one to which the monkeys had to then move the cursor. If the monkey achieved this and held the cursor on the new target for 300 ms, he received a juice reward.
Neuron to EMG joint cross-correlation.
After conducting preliminary analysis (Zhuang et al. 2010), we preferred JCCs to floating cross-correlation histograms (CCHs) (Moore et al. 1970; Schwartz and Adams 1995) because JCCs detected time-dependent changes in correlations with a 50-ms resolution, whereas CCHs required a 400-ms or larger floating window. Indeed, CCH represents correlation as a function of a single temporal parameter (lag). The data for these analyses are taken from a time window, which needs to be large enough to provide sufficient statistics. On the other hand, to capture dynamic changes, this time window is made relatively narrow (e.g., 350 ms in Schwartz and Adams 1995) and is slid along the time axis. Our examination showed that a sliding window of 400–1,000 ms would provide a sufficient statistical sample for our data, but the time resolution for detecting time dependency of correlations would be insufficient with such a large window. JCC resolves this problem by analyzing trial-by-trial correlation between neural signals sampled at particular intratrial times relative to task events. Furthermore, this method reveals certain features of cross-correlation that would be obscured by single-lag representations.
To compute JCCs, we analyzed self-timed touch pad press data from a total of two behavioral sessions in monkey A and used 1,000 behavioral trials per session. Two behavioral sessions of the center-out task in each of monkeys B and C were used, as well. For monkey B, we used 1,166 trials (145 ± 11 trials per reach direction) for one session and 605 (74 ± 6 trials per reach direction) for the other session. For monkey C, we used 1,211 trials (150 ± 7 trials per reach direction) for one session and 636 (80 ± 8 trials per reach direction) for the other session. Firing rates of individual cortical units were evaluated as time series of spike counts in 10-ms bins and then rescaled via z-scores. The resulting signal was low-pass filtered to below 20 Hz using a fourth-order Butterworth filter. EMG signals were also preprocessed. They were full-wave rectified, low-pass filtered with a fourth-order Butterworth filter with the same 20-Hz cutoff as for neuronal data and resampled to 10-ms sampling intervals. As a result, neuronal and EMG signals matched each other in sampling rate (100 Hz) and could be entered to the JCC algorithm.
Overall, our JCC method is very similar to the method known as JPSTH previously introduced for the analysis of correlation between pairs of neurons (Aertsen et al. 1989; Brody 1999; Brown et al. 2004; Cardoso de Oliveira et al. 2001; Gerstein and Perkel 1969), and a version of this analysis called the joint peri-event correlogram (Cardoso de Oliveira et al. 2001). We also considered trial-by-trial correlations between EMGs and neuronal rates as a function of their measurement times: EMG time, t1, and neuronal time, t2, both measured relative to a task event (Fig. 4A). To calculate JCC, we first selected a reference task event: the touch pad release onset in the timing experiment and peripheral target appearance in the arm reaching experiment. JCCs were calculated in a time frame around the task event, and the correlation was graphically represented as pixels (Fig. 4A). For the touch pad press task, the time frame we analyzed was 3.5 s before and 1.5 s after the release onset for both EMGs and neural units. For the center-out reaching task, the time frame analyzed was 0.5 s before and 1.5 s after time of peripheral target appearance on the screen. Average values of the rectified EMG and neural signals were calculated for 50-ms bins, with EMGs on the horizontal axis and neural units on the vertical axis of the JCC matrix. For each bin, Pearson's correlation coefficient was calculated to evaluate trial-by-trial cross-correlation between neuronal and EMG signals: (1)
where ρt1,t2 is the correlation coefficient between one unit and a single EMG channel at EMG lag t1 and neuronal lag t2. n(t2) is the vector (length equal to the number of trials) of neural activity at time lag t1 relative to a task event, and m(t1) is the vector of EMG activity at time lag t1 relative to the same task event for all trials. To construct the JCC matrix, this correlation calculation was done for each pair of neuronal and EMG time lags.
As shown in Eq. 1, each behavioral trial contributed a pair of numbers to every bin: readings of neuronal activity and EMG. The sample of these readings for all trials was then used to calculate the correlation coefficient. Thus JCCs evaluated the trial-by-trial correlation between the EMG and neuronal activity for each t1 and t2. These correlation data produced a two-dimensional (2D) JCC matrix. Similar matrices of correlation coefficients were previously used by Cardoso de Oliveira and colleagues (2001) to analyze functional connectivity between cortical neurons, but, to the best of our knowledge, this analysis has not been applied to cortico-EMG relationships.
Figure 4A shows a diagram explaining how to read a JCC graph. The diagonal (designated by a dashed line) represents correlation between simultaneous readings of the EMG and neuronal activity. Pixels above this line correspond to neuronal time bins that are later relative to EMG time bins. Conversely, pixels below the main rightward ascending diagonal indicate neuronal time bins that are earlier relative to EMG time bins. Figure 4A also shows a vertical band that corresponds to a fixed EMG time and different times of spike occurrences. This particular band was of interest to us because it corresponded to movement initiation time and an accompanying EMG burst.
As a statistical control, we entered randomized (Edgington 1995) EMG trials into the JCC algorithm (Figs. 4D and 5A, bottom). Correlation was calculated by using neural data from one trial and EMG data from a different trial. These EMG trials were shuffled without replacement (Fig. 4D). This randomization procedure, when repeated 400 times, provided a statistical sample for the estimation of confidence intervals. For each pixel in the JCC, we calculated the variance using randomized EMG-neural unit paired trials as a statistical sample. We used 95% confidence intervals to reject noise in JCCs based on this variance. Given this confidence interval, up to 5% of the JCC pixels could still represent values falsely classified as different from random. All JCC features described below were well above the chosen statistical threshold.
Comparison of JCC with STA.
The results obtained with JCC were compared with the measurement of correlation based on STAs. STAs were calculated from EMG data 0.2 s before to 0.2 s after each spike event (i.e., a 0.4-s interval); data points were sampled at 500 Hz. To analyze STAs during different time epochs relative to task events, spikes that occurred only during specific time epochs were used for calculation of the STA (Fig. 4B). Each of the four epochs was ∼0.5 s long: −1.25 s to −0.75 s, −0.75 s to −0.25 s, −0.25 s to 0.25 s, and 0.25 s to 0.75 s, inclusive, relative to touch pad release. We also calculated STAs with shuffled EMG trials (Fig. 4D), much like the shuffled JCCs (Fig. 5A, bottom). In this case, we calculated STAs with neural spikes during specific time epochs, but used averaged EMG activity in those time epochs from a different trial. This way, we obtained STAs that were corrected for randomized trials (Fig. 4C; red traces).
Note that the temporal resolution of JCCs for detecting changes in correlation was equal to JCC bin size (50 ms), whereas the resolution of STAs was limited by the epoch from which data were gathered (500 ms). These estimations of temporal resolution should not be confused with the capacity of STA to resolve the latency of rapid postsynaptic events, which can be 3–4 ms and is best detected when STA bins are small and the analysis window is long (Holmes and Spence 2004). Although such STAs are very good for detecting postsynaptic events and in this context can be described as having high temporal resolution, their temporal resolution for detecting time dependency of correlation may be poor when data are taken, regardless of task events from long time periods.
To specifically assess neural predictions of trial-by-trial variations of EMGs, we developed a linear decoder that predicted EMG amplitude on each trial from cortical signals. This trial-by-trial decoding should not be confused with continuous predictions of EMG from cortical activity (Morrow and Miller 2003; Santucci et al. 2005), where an overall correlation coefficient for continuous data is usually used as a measure of prediction accuracy. Trial-by-trial prediction accuracy for particular time points is rarely, if ever, considered.
To obtain EMG predictions from the firing of individual units, we applied a 2D linear regression, where EMG amplitude was predicted from firing rate for each JCC bin (see Fig. 9A, right): (2) (3)
where emg is the EMG signal during the time lag t1 across trials, is predicted EMG signal, A is the vector of regression weights of length N + 1 (N being total number of neural units) and n is a number of trials × (N + 1) matrix of neural activity for all neural units padded with ones to obtain parameters for offset. Equation 2 describes the prediction process used on test trials, whereas Eq. 3 describes how the regression weights are calculated using a training set of trials. In addition to applying Eqs. 2 and 3 to an ensemble of N neurons, we used the same equations to calculate single-neuron predictions (i.e., the case of N = 1).
Of 400 trials, 120 (30%) were used for training the model (Eq. 3). The parameters from the regression were then applied to the other 280 (70%) trials (Eq. 2). This process was repeated 100 times, each time with a random set of non-overlapping training and fitting trials. As a measure of prediction accuracy, we used a trial-by-trial correlation coefficient between the predicted and actual EMG amplitudes, which was calculated separately for each JCC bin. Figure 9A shows the average correlation coefficient for each JCC bin over the 100 fitting repetitions. Note that the prediction accuracies are presented in basically the same format as JCCs. Moreover, there is an excellent match between the values of JCC and prediction accuracies.
In a different prediction analysis, we used a 16-degree-of-freedom linear regression model sampled during a 0.75-s time window before release (see Fig. 9B) to predict the magnitude of the ensuing EMG burst for one neural unit at a time. The equation for training and fitting the linear regression model is the same as in Eqs. 2 and 3, except instead of nt2, we have tn: a number of trials × (number of time lags + 1) matrix of neural activity during the 0.75-s time window (consisting of fifteen 50-ms time lags) for a single neural unit, n. (4)
Here, M is the total number of trials used for fitting or testing and l to l + 15 are the 15 time lags during the time window of interest. A would now be a vector of length number of time lags + 1.
The regression was fit for each individual unit over a time of 0.75 s before the EMG burst. Four hundred trials were used: 120 for training and 280 for testing. Fitting and testing were repeated 100 times to obtain a mean correlation between actual and fitted EMG amplitudes for each neural unit. The histogram in Fig. 9C shows the distribution of prediction accuracies for different units.
In the self-timing experiment (Fig. 2), two sessions of simultaneous neuronal and EMG recordings were conducted in monkey A with 110 and 104 neural units (for each session, respectively) recorded in left (contralateral) M1, 94 and 84 in left PMd, 56 and 58 in right M1, and 96 and 85 units in right PMd. In the center-out reaching experiments (Fig. 3), a total of 19 M1 units were recorded simultaneously in monkey B in each of two recording sessions. In monkey C, 215 and 216 neural units were recorded simultaneously in each of two recording sessions (118 M1 and 97 S1 units in one session and 119 M1 and 97 S1 units in another session).
JCC analysis (Fig. 4A) revealed a variety of time-dependent correlation patterns (e.g., Figs. 5, 6, 9, and 10). These patterns can be seen in several representative examples for the self-timing task. Figure 2A (right) shows simultaneously recorded EMG traces of four muscles and activity of three M1 neural units contralateral to the working arm, recorded during the execution of the timing task. Peri-event time histograms (PETHs) for these neurons and average rectified EMGs for the same muscles are shown in Fig. 2B (left) aligned on touch pad release. Additionally, Fig. 2B (right) shows PETHs for the entire population of 356 cortical units recorded during one session (110, 94, 56 and 96 units recorded in contralateral M1, contralateral PMd, ipsilateral M1 and ipsilateral PMd, respectively). Clear EMG bursts occurred when the monkey released the touch pad (dotted lines), but there was also ongoing EMG activity during the delay period, especially evident in the wrist extensors. The extensors exhibited an EMG burst shortly before the touch pad release, followed by attenuated EMG activity, and another burst around the time when the hand returned to the touch pad. Wrist flexors exhibited bursts shortly before movement onset. These burst sequences in the forearm muscles are similar to the triphasic EMG pattern typical for ramp and hold movements, where a burst of activity in the agonist muscles (wrist extensors in our case) is followed by a pause in EMG activity in those muscles (Brown and Cooke 1990; Hannaford and Stark 1985; Marsden et al. 1983). The antagonist muscles (wrist flexors) then burst, and finally antagonists relax as the agonists reactivate. In our experiment, the amplitude of EMG bursts was clearly variable from trial to trial (Fig. 2A, right). Neuronal firing patterns also exhibited trial-to-trial variability (e.g., modulations around touch pad release in neurons illustrated in Fig. 2A).
To quantify cortico-EMG correlations for the timing task, we aligned JCCs on touch pad release (Figs. 5A and 6). The results were very similar for the alignment on touch pad press (not shown). Figure 5A shows JCCs for the M1 neuron marked as neuron “1” in Fig. 2A, right. JCCs of this neuron with wrist flexors and wrist extensors show several JCC features particularly well. These features can be described as bands of positive (red color) or negative (blue) correlations, each exposing a particular time-dependent pattern of cross-correlation. JCC bands have different orientations: notice diagonally and vertically oriented bands (marked in Fig. 5A).
The diagonal orientation of JCC bands corresponds to a classical cross-correlation pattern which is a function of the temporal lag between the two correlated signals, Δt = t1 − t2. (A diagonal line is described by the equation t1 − t2 = const.) In Fig. 5A, multiple diagonal bands of both positive and negative correlation are prominent around the time of touch pad release, t1, t2 = 0 s, and around t1, t2 = −3 s, corresponding to touch pad release that occurred on the preceding trial. The correlation signs are inverted for the flexors vs. extensors. Conspicuously, for both flexors and extensors, diagonal bands are present for negative Δt. Such negative lags correspond to neuronal modulations that occurred after changes in the EMG as if the neuron tracked the EMG modulations. Among the plausible explanations, these correlation peaks for which neuronal discharges lag the EMG may reflect cortical processing of movement-related reafferent signals (Barlow 1959; Brazier and Casby 1952; Perkel et al. 1967).
As explained in Fig. 4B, we converted JCC data into classical cross-correlations and calculated corresponding STAs and shuffle-corrected STAs for four task periods. The epoch 1.25 s to 0.75 s before touch pad release is indicated in blue, 0.75 s to 0.25 s before touch pad release is indicated in green, 0.25 s before to 0.25 s after touch pad release is indicated in yellow and 0.25 s to 0.75 s after touch pad release is indicated in orange. The results of this analysis for the neuron of Fig. 5A are shown in Fig. 5B. Notice that, in all cases, STAs for shuffled data (gray traces) were not flat, indicating the presence of comodulations of EMGs and neuronal rates that are not removed by the shuffling. The shuffle-corrected STAs (red traces) depict trial-by-trial comodulations, usually called functional or effective connectivity (Aertsen et al. 1989). It is clearly seen from Fig. 5B that JCC data and shuffle-corrected STAs are in good correspondence. However, JCCs of Fig. 5A are more useful than STAs for inspecting task-related changes in cross-correlations and detecting correlation features.
In Fig. 5A, vertical bands are prominent around t1 = 0 s, reflecting long-term correlations between neuronal activity and the EMG burst at the time of movement initiation. The presence of vertical bands indicates that the EMG burst amplitude was correlated with neuronal rate recorded as early as 1–2 s before the burst. Wrist flexors were negatively correlated with the neuron for neuron lag in the interval t2 = −2.0 to −0.2 s (blue vertical band). By contrast with this pattern, wrist extensors were positively correlated with the neuron for t2 during the interval −1.2 s to −0.2 s (red vertical band). Note that the vertical bands are only present in the JCCs of Fig. 5A around the time of touch pad release (t1, t2 = 0 s), and that they are completely smeared for the preceding touch pad release (t1 around −3 s), where there is no precise alignment on release onset.
To verify whether our JCC results were statistically sound, we conducted a shuffle test where single-neuron trials were matched to EMG trials in a randomized order (Fig. 5A, bottom) (Edgington 1995). After shuffling, correlation between the neuron and EMG traces was clearly low, and none of the original correlation features were present. After the pixels with statistically insignificant values of correlation were discarded, the JCC patterns reported here were mostly unaltered (not shown). Thus, unlike STAs, JCCs depicted only trial-by-trial correlations (i.e., functional connectivity) and did not reflect trial-independent comodulations.
To estimate the effect of shuffling across neuronal units, we characterized the presence of JCC features as JCC standard deviation. JCC standard deviation was calculated for all pixels in a JCC, and it characterized an overall variability. JCC standard deviation was then calculated after shuffling. The difference between the unshuffled and shuffled JCCs was highly significant across units: 0.039 ± 0.010, 0.041 ± 0.017, and 0.052 ± 0.022 (mean ± SD) before shuffling vs. 0.031 ± 0.002, 0.033 ± 0.011, and 0.038 ± 0.007 after shuffling in monkeys A (P < 0.001; Wilcoxon rank sum test), B (P < 0.001) and C (P < 0.001), respectively. Furthermore, the effect of shuffling was analyzed on a unit-by-unit basis. In this analysis, JCC standard deviation was compared with the standard deviation of shuffled JCCs using a randomization test (Edgington and Onghena 2007). Shuffling was performed 1,000 times to compute a statistical distribution and confidence intervals for the JCC standard deviation. We found statistically significant differences (P < 0.01) in 491/687 (71.5%), 38/38 (100.0%), and 424/431 (98.4%) neuronal units in monkeys A, B and C, respectively. To estimate the proportion of neural units exhibiting correlation bands depicted in Fig. 5A, we performed automated feature detection. Diagonal bands were classified as such if a 0.5-s-long band existed in the JCC that had a mean absolute value of correlation that was 2 SDs above that of the average 0.5-s-long diagonal bands. Similarly, vertical and horizontal bands were determined using the same 2 SD threshold over 0.5 s in either EMG time (for horizontal bands) or neuronal time (for vertical bands). Diagonal patterns were the prevalent feature among neural units. They were found for at least one muscle for 69.6% (149/214; 2 recording sessions) of contralateral M1 neural units and 50.0% (89/178) of contralateral PMd units. For the ipsilateral hemisphere, these frequencies were much lower: 28.1% (32/114) for M1 and 16.6% (30/181) for PMd. Vertical bands occurred in 36.0% (77/214), 20.2% (36/178), 21.1% (24/114) and 17.1% (31/181) of contralateral M1, contralateral PMd, ipsilateral M1 and ipsilateral PMd neural units, respectively. Horizontal bands near movement onset were observed for 10.7% (23/214), 3.4% (6/178), 0.0% (0/114) and 2.8% (5/181) of contralateral M1, contralateral PMd, ipsilateral M1 and ipsilateral PMd units, respectively.
For the center-out reaching task, diagonal bands were found for at least one muscle for 50% (20/40) M1 units recorded over two experimental sessions. Vertical bands were found for no units, and horizontal bands were found in 65% (26/40) of M1 units in monkey B. For monkey C, diagonal bands were observed in 77.6% (184/237; 2 recording sessions) of contralateral M1 units and 82.0% (159/194) of contralateral S1 units in monkey C. Vertical bands occurred in 46.0% (109/237) of M1 and 48.5% (94/194) of S1 units. Horizontal bands occurred in 29.5% (70/237) of M1 and 33.0% (64/194) of S1 units.
JCCs for antagonistic muscles.
Correlations, as revealed by JCCs for the self-timing task, did not necessarily have opposite signs for antagonistic muscles. For example, Fig. 5A shows that the neuronal activity was positively correlated with the EMGs of both wrist flexors and extensors throughout the delay period of the task (t1 and t2 in the interval −2 s to −0.2 s, indicated by dashed lines), as evident from the red-colored band adjacent to the diagonal. For both muscle groups, this band was broad (Δt on the order of 1 s) early in the delay and narrowed closer to movement onset. Notice that the narrow peaks occurred at different lags for these muscles (e.g., in Fig. 5B, correlation peaks for these muscles are separated by ∼100 ms for the peri-movement interval).
JCC patterns for two additional neurons are presented in Fig. 6. The JCCs shown in Fig. 6A (the same M1 neuron as “neuron 2” in Fig. 2A) are clearly different from the previous examples. This M1 neuron exhibited a negative correlation with wrist flexor EMGs, as can be seen from the band immediately adjacent to the JCC diagonal (Δt < 0.1 s). This negative correlation persisted during the delay period and strengthened around movement onset. During the delay period, a broad diagonal band indicated a longer-term negative correlation with the flexor EMG. The diagonal delay-period correlation with the extensor EMG was also negative, but weaker. Additionally, multiple patches of positive and negative correlations can be seen around touch pad release in JCCs for both wrist flexors and extensors, as well as for the biceps and triceps. The JCCs for the M1 neuron illustrated in Fig. 6B (the same as “neuron 3” in Fig. 2A) somewhat resemble the JCCs of Fig. 5A, particularly for the wrist flexors and extensors.
To analyze the correspondence between correlations of neurons to opposing muscles, we analyzed JCC segments for four task epochs (Fig. 4B). These segments were computed separately for JCCs of each neuron with the wrist flexors and the wrist extensors. Next, whether or not JCCs matched for the two muscles was determined using correlation (not to be confused with the other usage of correlation coefficient in this paper). This correlation coefficient evaluated the correlation between the JCC-derived cross-correlogram as shown in Fig. 4B. If the resulting correlation coefficient was significantly positive for a particular neuron (t-test; P < 0.01), it meant that this neuron has matching JCC pattern with both the wrist flexors and extensors. If the coefficient was significantly negative, it meant that the neuron reversed JCC patterns between the agonist and antagonist muscles. During the −1.25-s to −0.75-s epochs, 34/214 M1 neural units contralateral to the working arm (15.9%), 13/178 (7.3%) contralateral PMd units, 21/114 (18.4%) ipsilateral M1 units, and 44/181 (24.3%) ipsilateral PMd units displayed matching correlations. 17/314 (5.4%) contralateral M1 neural units, 4/178 (2.3%) contralateral PMd units, 4/114 (3.5%) ipsilateral M1 units and 15/181 (8.3%) ipsilateral PMd units displayed inverse correlations. During the −0.75-s to −0.25-s epoch, 48/214 (22.4%) contralateral M1 units, 43/178 (24.2%) contralateral PMd units, 31/114 (27.2%) ipsilateral M1 units, and 42/181 (23.2%) ipsilateral PMd units displayed matching correlations. 29/314 (9.2%) M1 neural units, 18/178 (10.1%) contralateral PMd units, 7/114 (6.1%) ipsilateral M1 units and 21/181 (11.6%) ipsilateral PMd units displayed inverse correlations. During the −0.25-s to 0.25-s epoch, 33/214 (15.4%) contralateral M1 units, 28/178 (15.7%) contralateral PMd units, 20/114 (17.5%) ipsilateral M1 units, and 36/181 (19.9%) ipsilateral PMd units displayed matching correlations. 51/214 (23.8%) contralateral M1 units, 27/178 (15.2%) contralateral PMd units, 24/114 (21.1%) ipsilateral M1 units and 29/181 (16.0%) ipsilateral PMd units displayed inverse correlations. During the 0.25-s to 0.75-s epoch, 12/214 (5.6%) M1 units, 26/178 (14.6%) contralateral PMd units, 22/114 (19.3%) ipsilateral M1 units, and 32/181 (17.7%) ipsilateral PMd units displayed matching correlations. 48/214 (22.4%) M1 units, 25/178 (14.0%) contralateral PMd units, 12/114 (10.5%) ipsilateral M1 units and 20/181 (11.0%) ipsilateral PMd units displayed inverse correlations.
JCCs and STAs for cortical populations.
To visualize JCC data for the entire ensemble of simultaneously recorded cortical units, we rendered a three-dimensional (3D) plot where 2D individual-unit JCCs were stacked into a 3D matrix with the dimensions represented by t1 (EMG lags relative to a task event), t2 (neuronal lags relative to the same task event) and number of neural units in the Z-axis (Fig. 7A). Each of the panels in Fig. 7C represents the average correlation over a time window along the main diagonal (blue plane of the lower panel in Fig. 7A). The time epochs labeled in this plot correspond to time relative to touch pad release onset (intersection of all three planes). It can be appreciated from this representation that practically all neural units in the recorded M1 population exhibited transient correlations with the wrist extensor EMGs.
The population plots of Fig. 7, B and C, show the evolution of lag-dependent STA (Fig. 7B) and cross-correlation (Fig. 7C, i.e., cross-correlation represented as a function of Δt) during the execution of the task. For each neural unit, lag-dependent cross-correlation at time t is represented by a slice through the JCC perpendicular to the main diagonal. For each of these epochs, a variety of correlation patterns characterized by the timing of positive and negative peaks can be seen. Notice that the timing of correlation peaks and troughs changes across the task epochs. Indeed, in Fig. 7, B and C, neural units are sorted according to the time of positive peak occurrence during touch pad release (from top to bottom, starting with the earliest peak). This orderly peak pattern (Fig. 7, B and C, third column) is distorted for the other task intervals.
For each JCC pattern shown in Fig. 7C, we calculated a matching shuffle-corrected STA. The degree of correspondence between the JCC and STA data was measured as correlation coefficient, RJCC-STA. The values of RJCC-STA for all contralateral M1 neural units during a single session are shown in Fig. 7D. Notice high RJCC-STA (i.e., matching STAs and JCC-derived curves) for many units. Across both sessions, RJCC-STA values were 0.164 ± 0.273 (mean ± SD) for −1.25 s to −0.75 s, 0.16 ± 0.34 for −0.75 s to −0.25 s, 0.29 ± 0.34 for −0.25 s to 0.25 s, and 0.18 ± 0.34 for 0.25 s to 0.75 s, all significantly different from zero (P < 0.001; Student's t-test) when averaged across all neural units and all muscles. These values were even higher when the top 50% neural units with high JCC modulations were selected: 0.21 ± 0.29 for −1.25 s to −0.75 s, 0.21 ± 0.36 for −0.75 s to −0.25 s, 0.43 ± 0.35 for −0.25 s to 0.25 s, and 0.26 ± 0.36 for 0.25 s to 0.75 s (P < 0.001). In other words, for the units with high cross-correlation with EMG, STAs and JCC-derived curves matched well. This trend was true when averaged across brain areas: areas with lower cross-correlations (i.e., ipsilateral M1 and PMd), displayed lower RJCC-STA. When averaged over all task epochs, RJCC-STA for the different brain areas was 0.30 ± 0.52 for contralateral M1 neural units, 0.22 ± 0.51 for contralateral PMd units, 0.09 ± 0.50 for ipsilateral M1 units, and 0.14 ± 0.47 for ipsilateral PMd units (P < 0.001).
Changes in correlation width.
It can be noticed in Figs. 5–7 that the correlation curves have more narrow peaks (or troughs for negative correlation) during touch pad release (third column) compared with the other periods. Figure 8 shows a quantitative analysis of this effect. The width of the correlation curve, defined as the time over which correlation stays above one-half of the midpoint between correlation peak and trough (full width at one-half of maximum, FWHM), was on average the shortest during the peri-movement period (−0.25 s to 0.25 s) and the longest for the early delay period (−1.25 s to −0.75 s). This effect was seen for forearm (flexors and extensors) and arm (biceps and triceps) muscles and for all cortical areas (contra- and ipsilateral M1 and PMd) (P < 0.01, Kruskal-Wallis test with α-level of 0.01 post hoc comparison). Correlation width was thus narrower during movement when EMG activity was modulated more synchronously with neurons.
Vertical band and single-trial predictions.
The origin of the long-term cortico-EMG correlations revealed by the vertical bands (e.g., Figs. 5 and 6) was not immediately clear from the JCC analysis alone. Figure 9 depicts an extended analysis for this particular JCC feature. Figure 9A shows two neurons with vertical JCC bands. One neuron exhibited positive vertical-band correlation for the data points when neuronal activity preceded the EMG of wrist flexors and negative correlation when neuronal activity succeeded the EMG (Fig. 9A, top left). Another neuron exhibited multiple vertical bands in the JCC for the wrist flexors, which also reflected positive or negative correlation, depending on the pattern location within the JCC (Fig. 9A, bottom, left).
To clarify the origin of vertical bands, such as those in Fig. 9A, we split trials into a subset with low-amplitude EMG bursts and one with high-amplitude EMG bursts. First, the time bin with the highest EMG amplitude was found for each muscle group across trials. At this time bin, we categorized trials during which the EMG burst amplitude exceeded that of the average burst amplitude across trials as “high-EMG burst” trials. Likewise, trials with burst amplitudes lower than the average were categorized as “low-EMG burst” trials. Average rectified EMGs (Fig. 9B, top) and PETHs (see PETH for an example neuron Fig. 9B, bottom) were calculated separately for these types of trials. This analysis revealed that low-EMG burst trials were accompanied by increases in EMG activity during the delay period of the task, whereas there were no EMG modulations during this period in high-EMG burst trials. From Fig. 9B (bottom) one can also see that neuronal activity during the delay period was also different for these types of trials. The neuron illustrated in Fig. 9B (bottom) was more active during low-EMG burst trials than during high-EMG burst trials, hence a negative correlation for the vertical band (Fig. 9A, bottom, left). Thus the occurrence of JCC vertical bands was explainable by the monkey's motor behavior. In a portion of trials, the monkey started to move the hand during the delay interval while holding the hand on the touch pad. In these trials, less additional muscular force was needed to release the touch pad, so the EMG burst was lower in amplitude. Motor cortical neural units tracked this motor behavior, which resulted in a long-range correlation between the units and the EMG.
Furthermore, we analyzed single-bin predictions of EMGs for each individual unit (Fig. 9A, right). To calculate these predictions, we applied Eqs. 2 and 3 to a single neuron (N = 1). Accuracy of these predictions, evaluated as correlation coefficient, was in excellent correspondence with JCC values (compare left and right panels in Fig. 9A). Indeed, the presence of positive or negative correlation for a particular JCC bin translated into a high accuracy of EMG predictions for that particular bin.
A particular example of such good correspondence between JCCs and EMG predictions is shown in Fig. 9C. Here, the population of M1 neural units was split into the subpopulation with vertical JCC bands (i.e., neuronal activity positively or negatively correlated with EMG bursts) and the subpopulation without vertical bands (i.e., no correlation). When trial-by-trial predictions of EMG burst amplitude were calculated separately for these subpopulations, the accuracy of individual-unit predictions was clearly higher for the subpopulation with vertical bands (mean ± SD for the R of individual-unit prediction of 0.053 ± 0.008 for units with vertical bands, 0.033 ± 0.007 for units without vertical bands, P value ≪ 0.01). While these R values are low, we note that the goal of this analysis was not to produce practical BMI predictions, but rather to use individual-unit predictions as a measure of neuronal encoding.
It should be also noted that our ability to predict EMGs from cortical activity was not necessarily a result of a causal relationship. The nature of correlation underlying this prediction could be quite complex. For instance, as shown in Fig. 9D, EMGs of different muscles were correlated themselves. This figure shows JCCs between pairs of EMG channels, and it can be seen that the time-varying effects seen between neurons and muscles are evident between muscles, as well.
Overall, we observed that particular JCC features, such as vertical bands, could be used to detect motor behaviors and associated cortico-EMG relationships, which would otherwise go unnoticed. Furthermore, the presence of positive or negative correlation for particular JCC components indicated that EMGs could be predicted from neuronal activity for the corresponding values of t1 and t2.
Reaching-task JCCs and directional tuning.
For the center-out reaching task (Fig. 3), JCC analysis revealed not only time-dependent correlation patterns between cortical units and arm EMGs, but also time-dependent directional properties of those correlations.
Figure 3 illustrates example EMG and neuronal recordings (Fig. 3B), population PETHs for different directions of movement (Fig. 3, C and D, left) and average EMG traces for different directions (Fig. 3, C and D, right). Data from both monkey B and monkey C are shown. Movement-related neuronal modulations and EMGs are clearly directionally tuned.
Figure 10 shows a novel directional-tuning plot. In this analysis, JCCs are constructed for eight directions of reaching movements. The reach direction for any particular trial was picked from a random continuous distribution across all 360°. For our analysis, however, we grouped reach directions into one of eight ranges separated at 0°, 45°, 90°, 135°, 180°, 225°, 270° and 315°. The times for the EMG channel, t1, and the neural unit, t2, are measured with respect to peripheral target appearance. The JCCs show that four representative M1 neurons (three from monkey C and one from monkey B) were correlated with the wrist flexor EMG in monkey B and the triceps EMG in monkey C, both before and during arm reaching movements. The JCC tuning patterns are different for each illustrated neuron. Thus the neuron shown in Fig. 10A had the highest correlation with the EMG for the 315° target location. For 135°, the opposite direction, the correlation peaked with a negative sign (i.e., blue color in the JCC). For the neuron illustrated in Fig. 10B, the JCC tuning pattern was bimodal, with peak positive correlations at 45° and 180°. The JCC peak for the neuron shown in Fig. 10C was prominently negative at the 0° target location. Finally, the JCC for the neuron of Fig. 10D, recorded in monkey B's M1, has peak positive tuning at 225°.
In addition to JCCs for different movement directions, Fig. 10 shows polar plots for directional tuning of each neuron, EMG and JCC calculated for the period 0.5 s before to 2.5 s after target onset. Notice that the tuning curves for different signals and their correlations do not match each other. For example, the preferred direction for the neuron in Fig. 10C is at 330°, the EMG's preferred direction is at 250°, and the negative correlation of JCC peaked at 0° (positively toward 180°), the direction for which minimal modulations in both neuronal rate and EMG occurred.
Moreover, directional tuning curves for each parameter exhibited time dependency. This is clear from the diagrams at the bottom of each figure panel of Fig. 10, which show the mean vector for the neuron, EMG and JCC as a function of time relative to target onset. It can be seen that both the amplitude and the direction of the mean vectors change with time.
JCCs directional tuning for cortical populations.
Figure 11 illustrates the diversity of directional tuning of cortico-EMG correlations for all neural units recorded in monkeys B and C. Figure 11A shows color-coded directional tuning curves calculated for the period from 0 s to 0.5 s after target onset. In these plots, rows represent JCC tuning curves for different neural units, and different blocks correspond to muscle groups. Overall, a diversity of preferred directions was represented by different neuron-muscle pairs. This fact is illustrated in Fig. 11B, which shows the corresponding distributions of mean vectors as polar plots. To analyze the time dependency of JCC tuning for the entire neuronal population, we plotted their average correlation to each of the three or four muscles at each reach direction and over time (Fig. 11C). It can be seen that the neural population is more tuned to the wrist extensor and triceps muscle groups than the wrist flexors for all reach directions in monkey C, and that units are on average the least tuned to the triceps in monkey B. Neurons are more correlated to the wrist extensors and especially the triceps in monkey C during and after target onset.
JCCs directional tuning for cortical populations.
While directional tuning in terms of neural units, EMGs and JCCs were variable, they were not completely independent. More experiments examining different muscles and different cortical locations will be needed to clarify this relationship in detail. As a representative example from our data, Fig. 12 presents an analysis of directional tuning in monkey C. Directional tuning curves were calculated based on neuronal rates and JCCs of the same neural units with all four muscle groups recorded for the interval 0–0.5 s relative to target onset. The difference between the preferred direction for the neural unit and for the JCC was evaluated as the absolute value of the angle between the corresponding mean vectors. In this figure, each dot represents the difference between neural unit and JCC tuning direction plotted against the tuning direction for each unit. The red line indicates the tuning direction of the EMG signal for that particular muscle. In some pairs, there are obvious clusters in the data, particularly in the panel analyzing triceps vs. S1 unit correlations. Here, S1 units that closely match the EMG correlation direction display much more correspondence between their tuning and the JCC tuning directions. The widely varying tuning directions of M1 units also show that congruent tuning of units and muscles does not necessarily entail a congruent tuning direction of the JCC. Results in monkey B were less conclusive due to the lower number of units recorded (not shown). Still, clustering appeared in neural unit vs. JCC tuning with the wrist flexors. Here, differences in neural unit and JCC tuning were lower for units with tuning similar to the EMG tuning direction. Neural units in monkey B also showed sparse representation around 75°.
Correlation width for task epochs and directions.
As with the touch pad press task, the width of cross-correlation peaks was not static over the course of each trial during a center-out task. In monkey C, we found that the FWHM of the correlation peaks were narrowest during the time window 0–0.5 s after target appearance when averaged across all neural units, for all reach directions, considering all muscles. Mean of FWHM was 0.55 ± 0.37, 0.51 ± 0.43, 0.46 ± 0.33 and 0.53 ± 0.37 during the task epochs −1 s to −0.5 s, −0.5 s to 0 s, 0 s to 0.5 s and 0.5 s to 1 s with respect to target appearance, respectively. These differences were significant (P < 0.01) between the 0-s to 0.5-s task epoch and the other three epochs, but not significantly different between the other three task epochs. We also found that these dynamic changes in correlation width were directionally tuned. The 0-s to 0.5-s task epoch was only significantly different (P < 0.01) from the 0.5-s to 1-s task epoch in the 45°, 0°, and −45° reach directions, although it was significantly different from both pretarget appearance task epochs for all reach directions. In addition, overall correlation width was significantly narrower for the 0 and −45° reach directions across task epochs than the 225° and 270° reach directions. The 45° direction showed significantly narrower correlation widths than at 270°, even though the correlation width was not significantly different between any other pairs of reach directions.
In this study, we examined time-dependent changes in the correlation between cortical neuronal activity and EMGs of arm muscles in rhesus monkeys, while these animals performed motor tasks with their hands and arms. Cortical neuronal activity was recorded from chronically implanted multielectrode arrays placed in M1, PMd and S1. JCCs were employed to examine time-dependent changes in the correlation patterns between neuronal firing and EMGs of the forearm and arm muscles. This analysis revealed widespread correlations between cortical ensembles and multiple muscles. Practically every recorded cortical neural unit was correlated with EMGs of multiple muscles at some moment in time. Moreover, the JCC analysis clearly revealed that cortico-EMG relationships changed significantly with the passage of task epochs. Generally, these correlations could not be described only as a function of the lag between neuronal discharges and EMGs, but rather reflected more complex relationships, such as dependency of peri-movement EMG bursts on the delay-period activity of the neural units (vertical bands). When monkeys made center-out reaching movements with their arms, JCC patterns were directionally selective. Overall, these findings reveal a highly dynamic correspondence between cortical ensembles and spinal motoneurons, some of which was not apparent before we employed the JCC method introduced here.
These results add novel insights to previous studies of neuronal-EMG correlations, an area of research that started several decades ago. In the 1950s, correlation techniques were applied to electroencephalograms to reveal common activity in disparate regions of the brain as well as encoding of sensory inputs (Barlow 1959; Brazier and Casby 1952). In the 1960s, with the advent of single-unit recordings in behaving animals, correlation analysis was applied in an attempt to detect connectivity between neurons (Perkel et al. 1967). Similar approaches continue to be used to this day (Aertsen et al. 1989; Nevado et al. 2012). The now classic STA approach, which effectively represents cross-correlation between single neurons and EMGs, was introduced to assess functional connectivity of spinal motoneurons (Fetz et al. 1989; Jackson et al. 2006; Mantel and Lemon 1987; McKiernan et al. 2000; Miller et al. 1993; Morrow and Miller 2003; Schieber and Rivlis 2007). This and similar methods have been used to analyze the inputs to the spinal cord from the cortex (Carmena et al. 2003; Davidson et al. 2007; Fuglevand 2011; McKiernan et al. 2000; Morrow et al. 2007; Schaefer et al. 2009; Schieber 2011), red nucleus (Miller et al. 1993) and reticular formation (Stuphorn et al. 1999).
The analytic approach employed in the present study, JCC, is very similar to the JPSTH method introduced in the late 1960s for the analysis of time-dependent correlations between neurons (Aertsen et al. 1989; Cardoso de Oliveira et al. 2001; Gerstein and Perkel 1969; Vaadia et al. 1995). However, as far as we can tell, no studies have yet applied the technique of JPSTH to cortico-EMG correlations, despite the fact that these correlations are known to be time dependent, as previously shown using a floating window method (Schwartz and Adams 1995).
JCCs and STAs.
The results obtained with JCCs were generally consistent with STAs derived from the same data. However, STA analysis encountered several problems when applied to the same data. The first of these problems became evident when we calculated STAs for randomly shuffled EMG trials. Although the EMG and neuronal records were mismatched after this shuffling, their STAs were nonflat, indicating a significant proportion of comodulations related to common task-related activity in addition to individual trial-dependent correlation and commonly called effective connectivity (Aertsen et al. 1989). Common task-related activity may occur even if a neuron and a muscle are not directly interconnected, but each responds to task events. After we corrected STAs by subtracting the STAs for shuffled trials, we obtained a curve that reflected an extra correlation on top of such common modulations. The results obtained with these shuffled STAs corresponded well to the JCC data. The JCCs did not contain any common-modulation components, because they were based on trial-by-trial correlations.
The second problem related to using STAs to analyze dynamic changes in cross-correlation was the requirement to use a relatively wide data window [500 ms in our study; 350 ms in the analysis of cross-correlation of Schwartz and Adams (1995)]. This window was large compared with our 50-ms JCC bins. Therefore, JCCs employed better resolution for analyzing the fine time dependency of cross-correlations between cortical units and EMG. Still, 50-ms bins were too large for analyses of short-latency postsynaptic effects, usually conducted using STAs. For example, STAs resolve postspike activations and suppression with a latency of 6–8 ms (Holmes and Spence 2004). In this study, we did not reduce the size of JCC bins for two reasons. First, we did not have enough behavioral trials to gather sufficient statistics. Second, we used surface EMG recordings, whereas intramuscular EMG recordings are more suitable for revealing short-latency postspike effects (Holmes and Spence 2004; Schieber and Rivlis 2005). Having noted this, there is no principal limitation for making JCC bins smaller, and short-latency effects can be studied with this method in the future.
The third issue is a general convenience of using STAs for time-dependency analysis. With JCCs, task-related correlation patterns are readily visible in the plot. With STAs, an appropriate window size must first be selected, and then this window must be moved along the data to detect changes in cross-correlation. If these parameters are adjusted correctly, and shuffle-corrected STA frames are stacked together and displayed as a color plot, the resultant graph would look very similar to the JCC for the same data.
Time dependency of cross-correlation.
We found JCCs to be very effective for evaluating time-dependent correlations between recorded activity of cortical units and EMG activity of the arm and forearm muscles. Using this approach, we observed clear time-dependent correlation changes in two distinct experimental paradigms. JCCs visualized more data compared with the traditional approach that expresses the correlation between a neuron and an EMG as a function of the temporal lag between the two signals. The lag for which the correlation reaches a peak value is presumed to correspond to the signal transmission time (also called latency) from the neuron to the EMG. The existence of a correlation peak, often called postspike facilitation (PSF), or trough, called postspike suppression (PSS), at a latency of less than 9 ms is considered evidence that the neuron projects to the motoneuronal pool of the examined muscle through a monosynaptic or oligosynaptic connection (Schieber and Rivlis 2005). PSF is interpreted as being associated with excitatory connections (Buys et al. 1986; Fetz and Cheney 1980; Schieber and Rivlis 2005), whereas PSS is interpreted as an inhibitory connection (Kasser and Cheney 1985). The muscle for which the correlation effects are the strongest is called the target muscle for a given neuron. This STA analysis is equivalent to cross-correlation analysis (Lemon and Mantel 1989), but, as explained above, it reflects both common modulations and effective connectivity. The data for STA analysis are usually taken from long, continuous recordings, even if there are clear nonstationarities in the recording data, such as those related to the performance in a motor task (Griffin et al. 2008). As such, STAs often exhibit features different from short-latency peaks or troughs (Poliakov and Schieber 1998), many of which probably reflect common modulations of both the neuron and the EMG signal.
Our present JCC analysis did not resolve short-latency effects (because of 50-ms bins), but easily revealed dynamic changes in effective connectivity between cortical neurons and EMGs. In this context, we obtained clear evidence that there was no fixed correlation relationship for practically any cortical unit-EMG pair. Rather, such correlations changed over time with respect to a particular task event. Several previous studies came to similar conclusions. McKiernan et al. (2000) reported only a weak relationship between short-latency postspike effects (i.e., a relatively stable relationship based on direct connectivity) and more broad neuron-muscle covariations. Another study (Schwartz and Adams 1995) reported time-dependent changes in cortico-EMG cross-correlation during a sinusoidal tracking task. Zhou and Fuster (2000) reported changes in STAs after changes were made in the behavioral task. Morrow et al. (2007) concluded in their study that the relationship between neural and EMG activity was stable over time, but their time span of analysis was on the order of minutes and not gated to task events. None of these previous papers used an approach similar to our JCC method, so it would be of interest to apply JCCs to their data.
Types of JCC patterns.
In many cases, we observed lag-dependent cross-correlation patterns depicted by diagonal JCC features (several 50-ms bins on both sides of the major diagonal). In the literature, such lag-dependency is often interpreted in terms of relatively simple circuits in which cortical neurons maintain a fixed functional relationship with α-motoneurons through a monosynaptic or polysynaptic projection (Fetz and Cheney 1987). Additionally, lag-dependent correlations may occur through indirect effects. For example, this effect could emerge if a given cortical motoneuron does not innervate the muscle being recorded, but instead innervates a different muscle that acts in synchrony with the former muscle (McKiernan et al. 2000). It is also possible that some of the lag-dependent correlations observed in our study reflected the cortical unit's responses to proprioceptive feedback generated by muscle contractions. Indeed, EMG patterns in the forearm and arm muscles were consistent with the classical triphasic pattern of activation in agonist and antagonist muscles. The magnitude of EMG bursts in the agonist has been shown to be correlated with the antagonist's burst magnitude in normal subjects, but not in patients deprived of proprioceptive and cutaneous feedback (Forget and Lamarre 1987), indicating that the nervous system uses sensory feedback to modulate sequential motor bursts.
In addition to “classical” patterns, we also observed JCC features suggesting that correlations between neural units and EMG signals could be significant across a time period as long as 1 s. These unique patterns, such as vertical bands in the JCC matrix, were dependent on data alignment on task events (e.g., touch pad release). Imprecise alignment of the neural unit and EMG signals with such task events would result in smearing of such patterns.
Most generically, JCC features can be described as patches of negative and positive correlation. We classified these patches as diagonal, vertical and horizontal bands. For diagonal bands, we systematically observed that correlation peaks were narrower in the temporal vicinity of movement onset and broader at greater time lags. These changes in peak width were likely related to the level of synchrony in neuronal population activity. Indeed, many neural units fired in synchrony around movement onset. With such widespread synchrony, many units worked together as part of a circuit to produce the observed EMG bursts. Since cortical units could be correlated with each other, the resulting JCC could exhibit stronger and temporally compact correlation between a given unit and the EMG signals. Furthermore, broader correlation peaks in the absence of movement likely resulted from slower, weaker and less synchronous modulation that engaged cortical and spinal neurons during posture maintenance. Additionally, the level of neuronal activity probably reflected the expression of correlation peaks (Dietz 2003). Cohen and Kohn (2011) noted that subthreshold membrane potentials and their correlations could be masked by spiking threshold. This effect results in the dependency of the correlation on the firing rate: when the mean membrane potential is far below threshold, responses are weak, and many of the shared membrane potential fluctuations are unobservable in the spiking responses.
Directional tuning of cross-correlation.
In the center-out reaching task, we found that correlation patterns between pairs of neural units and EMGs were clearly dependent on the direction of reach. Interestingly, the preferred direction obtained for the neural unit-EMG correlation was generally unrelated to the preferred direction exhibited by the cortical unit used in the analysis. This finding bears similarity with the observations of Griffin et al. (2008), who found substantial disparities between the magnitude of PSF in neural unit-muscle pairs and the degree to which neuronal and EMG modulations matched each other in a reaching task. Griffin et al. showed that the peaks in cortico-motoneuronal cell activity very rarely coincided with peaks in EMG activity in the muscles those cells influenced.
Given a large number of factors that affect neuronal activity and EMGs and the linkage between them (i.e., the fact that the connectivity between a cortical neuron and spinal motoneurons is superimposed onto projections from other brain regions), it is not surprising to find nonmatching directional tuning in terms of firing rate and in terms of cross-correlation, as in our present study, as well as nonmatching occurrences of neuronal and EMG modulations as in Griffin et al. (2008).
Davidson et al. (2007) and Schieber (2011) called a similar disparity between activity of M1 neurons and EMGs, “dissociating motor cortex from the motor.” In their experiments, monkeys were rewarded for synchronous activation of specific motor neurons and arbitrary muscles. Monkeys were able to achieve the required comodulation by rapidly modifying functional connectivity of cortical neurons with α-motoneuron pools. Previously, the capacity of cortical neurons to disengage their modulations from the activity of target muscles was shown by Fetz and Finocchio (1971).
Time-dependent changes in neural unit-EMG correlations and their variability across units were especially striking in the population diagrams. The width, timing and even the sign of correlation peaks changed dramatically along the period of a single task trial. These findings indicate that the correspondence between cortical populations and spinal motoneurons is distributed, task and context dependent and highly dynamic.
A critic might argue that with 50-ms binned JCCs; we did not detect true connectivity between cortical populations and spinal motoneurons. For example, a neuron could become correlated with a given muscle indirectly if another muscle is correlated with the one to which the neuron is directly connected. Note, however, that JCC tends to attenuate indirect connectivity effects because it reflects effective connectivity rather than task-related comodulation. Furthermore, this issue could be tackled in the future by simply reducing the bin size used to calculate the JCCs.
Cortical unit-EMG correlations were observed for both M1 and S1 neurons, often with similar JCC patterns. This suggests that representation of muscle activity can be obtained from both sensory and motor cortical areas. The rich interconnections between M1 and S1 provide a neural circuitry landscape for such distributed representation (Pandya and Kuypers 1969). It has been noted since the mid-1900s that somatosensory and motor areas often exhibit similar physiological properties (Lilly 1956). More recent studies utilizing multielectrode recordings (Carmena et al. 2003; Schwarz et al. 2014), optical imaging and two-photon microscopy (Petreanu et al. 2009) have provided further evidence for corresponding activity patterns in these areas.
Previously, we argued that the representation of information by cortical neuronal ensembles in both M1 and S1 is highly distributed and dynamic (Nicolelis and Lebedev 2009). Here we observed a similar highly dynamic relationship between cortical neuronal ensembles and spinal motoneurons, whose output activity was recorded indirectly via surface EMG sampling. On a conservative note, the exact connectivity between the cortical spinal networks cannot be elucidated from JCC analysis alone. JCC patterns may reflect direct connectivity, common inputs or feedback loops. Each of these possibilities can be addressed in future studies using appropriate experimental designs.
Without making strong statements about the underlying connectivity, we note that several principles of neuronal ensemble physiology proposed by our laboratory (Nicolelis and Lebedev 2009) hold true for cortico-EMG relationships. For example, according to the distributed coding principle, any behavioral parameter is represented by multiple brain areas. This idea has also been investigated by Rathelot and Strick (2006). In their study, EMGs of arm muscles were correlated with neuronal ensembles in several cortical areas, such as M1 and S1. In another study, they found that different arm muscles had overlapping maps in M1, and most cortico-motoneuronal cells innervating these muscles were located in the anterior bank of the central sulcus (Grillner et al. 1981). We did observe higher correlations for the areas contralateral to the working arm, but significant correlations were also found for the ipsilateral hemisphere. Moreover, in the contralateral hemisphere, we did not find subpopulations of M1 or S1 neural units specialized for interacting only with particular muscles. Instead, the same ensemble represented several muscles, which is consistent with previous studies that employed STA (Holmes and Spence 2004; Lemon et al. 1986; Lemon and Mantel 1989). In addition, the single neuron insufficiency principle states that a single neuron is limited in encoding a given behavioral parameter. This principle is clearly valid for cortico-EMG relationships, because we did not observe that single neurons could mirror the activity of a particular muscle. The correlation between neural units and EMGs is strongly dependent on task interval, which in BMI applications makes it impossible to decode EMG from the activity of a single unit.
The multitasking principle proposed by us previously (Nicolelis and Lebedev 2009) indicates that a single neuron's firing can contribute to the encoding of several behavioral parameters. In the present study, we observed that firing of a single neuron was correlated simultaneously with EMGs obtained from several muscles. The applicability to the cortico-EMG relationships observed here and other principles of neural ensemble physiology proposed by us (degeneracy, plasticity, conservation of firing and context) would be an interesting subject for future studies.
Although the results of this study are insufficient to strongly claim that EMGs are generated by a highly distributed cortico-spinal network (rather than, for example, labeled lines), we prefer this interpretation.
JCC methods and BMIs.
The field of BMIs is expanding rapidly, and a variety of methods for decoding cortical signals have been proposed. One method of obtaining a control signal for BMIs is to decode EMG signals from cortical recordings, and several groups are actively exploring this decoding paradigm (Fagg et al. 2007; Nazarpour et al. 2012; Pohlmeyer et al. 2007; Santucci et al. 2005). One advantage EMG decoding has over kinematic decoding is that EMG signals are related to movement dynamics. Movement dynamics take into account how the actuator interacts with the environment, which is a key milestone for BMI systems that aim to reach clinical relevance.
These BMI studies often describe cortical neurons as encoders of EMG patterns, and, as such, characteristic correlations between neurons and EMG signals would provide the type of signal that is desirable for BMI control of an end effector. These correlations do not necessarily need to imply a true neuron-to-muscle connection as long as correlations persist throughout the behavior. Correlations revealed by JCCs can potentially reveal new neuron-EMG relationships that could be utilized in future BMI designs. In particular, our data showed that cortico-EMG correlations varied significantly as a function of time. Therefore, potential BMI applications will have to take this time dependency into account while designing robust neuroprosthetic applications.
In conclusion, we propose that JCCs may be useful to visualize the time-varying relationship between neurons and EMGs (as well other effectors) to characterize the performance of BMI decoders and plastic changes during long-term BMI control. As we know that correlations between neurons and EMG/effector activity change, depending on time relative to certain task events, it should be possible to use this information to develop more accurate decoding algorithms. For example, the decoder can take into account the effector's state in time to choose which neurons should be weighted more heavily at each temporal interval. In addition, neural plasticity may occur during BMI performance, such that the same set of neurons may contribute differently to the motor output over a given session. Without continuous updates to the decoder, the information present in neural recordings may be processed suboptimally during BMI control. Here, the JCC analysis was able to reveal temporal changes in neuronal-EMG correlations. This suggests that our analytic tool could be used to characterize neural plasticity in a quantitative way. As such, real-time JCC analysis could be employed to explore a sample of trials obtained by a BMI operation to continuously update the contributions of each individual neuron (and the entire neuronal ensemble recorded) to the task at hand. The obtained information could then be used to update the BMI decoding model to maximize its accuracy.
Research reported in this publication was supported by the National Institute of Neurological Disorders and Stroke of the National Institutes of Health (NIH) under award no. R01-NS-073952 and the National Institute of Mental Health of the NIH under award no. DP1-MH-099903 to M. A. L. Nicolelis.
The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.
No conflicts of interest, financial or otherwise, are declared by the author(s).
Author contributions: K.Z.Z., M.A.L., and M.A.L.N. conception and design of research; K.Z.Z. and M.A.L. performed experiments; K.Z.Z., M.A.L., and M.A.L.N. analyzed data; K.Z.Z., M.A.L., and M.A.L.N. interpreted results of experiments; K.Z.Z. prepared figures; K.Z.Z., M.A.L., and M.A.L.N. drafted manuscript; K.Z.Z., M.A.L., and M.A.L.N. edited and revised manuscript; M.A.L.N. approved final version of manuscript.
We are grateful to Dragan Dimitrov and Laura Oliveira for conducting implantation surgeries; Gary Lehew and Jim Meloy for engineering the experimental setup and the multielectrode arrays; Wenwen Bai, Joseph O'Doherty, Tamara Phillips and Andrew Tate for experimental support; Susan Halkiotis for administrative support and editing the manuscript.
- Copyright © 2014 the American Physiological Society