## Abstract

Interactions among groups of neurons in primary motor cortex (MI) may convey information about motor behavior. We investigated the information carried by interactions in MI of macaque monkeys using a novel multielectrode array to record simultaneously from 12–16 neurons during an arm-reaching task. Pairs of simultaneously recorded cells revealed significant correlations in their trial-to-trial firing rate variation when estimated over broad (600 msec) time intervals. This covariation was only weakly related to the preferred directions of the individual MI neurons estimated from the firing rate and did not vary significantly with interelectrode distance. Most significantly, in a portion of cell pairs, correlation strength varied with the direction of the arm movement. We evaluated to what extent correlated activity provided additional information about movement direction beyond that available in single neuron firing rate. A multivariate statistical model successfully classified direction from single trials of neural data. However, classification was consistently better when correlations were incorporated into the model as compared to one in which neurons were treated as independent encoders. Information-theoretic analysis demonstrated that interactions caused by correlated activity carry additional information about movement direction beyond that based on the firing rates of independently acting neurons. These results also show that cortical representations incorporating higher order features of population activity would be richer than codes based solely on firing rate, if such information can exploited by the nervous system.

A major problem in understanding higher brain function is determining how complex neuronal representations of sensory and motor information are constructed from the activities of individual neurons. Higher order sensory and motor representations appear to emerge from the firing of neuronal assemblies, but it has yet to be determined whether spatial and temporal interactions contribute to these representations. Interactions clearly occur within the widespread network of synaptically modifiable horizontal connections in which cortical neurons are embedded (Huntley and Jones, 1991; Hess and Donoghue, 1994; Rioult-Pedotti et al., 1998). Although controversial, various forms of interactions among cortical neurons, including firing synchrony and local oscillations, have been considered as key to forming higher level representations (Singer, 1995; Vaadia et al., 1995; DeCharms and Merzenich, 1996; Riehle et al., 1997; Hatsopoulos et al., 1998b).

The representation of information by populations of neurons has been extensively studied in both sensory (Nicolelis et al., 1995; Meister, 1996) and motor areas of cortex (Georgopoulos et al., 1988;Schwartz, 1993; Tanaka, 1994). A common observation is that the firing rates of single cortical neurons during repeated presentations of the same stimulus vary considerably (Henry et al., 1973; Tomko and Crapper, 1974; Heggelund and Albus, 1978; Rose, 1979; Dean, 1981; Lee et al., 1998; Shadlen and Newsome, 1998). Whether this variability is meaningless noise, reflecting the stochastic nature of neurons, or carries information is widely debated (Richmond and Optican, 1990; Bair and Koch, 1996; Buracas et al., 1998; Lee et al., 1998; Oram et al., 1998; Shadlen and Newsome, 1998; Zhang et al., 1998).

The firing rate of primary motor cortex (MI) neurons broadly reflects the direction of arm-reaching movements (Georgopoulos et al., 1986). The broad tuning of individual cells to movement direction coupled with their firing rate variability (Lee et al., 1998) have inspired theories of directional coding in which an accurate representation emerges as a result of averaging over a large number of cells (Georgopoulos et al., 1992; Salinas and Abbott, 1994; Sanger, 1994; Tanaka, 1994). Although these methods are quite successful in deriving movement direction, response variability was treated as a simple random noise process. This assumption is called into question both because individual MI neurons nearly always represent additional motor variables (Kalaska, 1988; Fu et al., 1993) and because weak correlations have been encountered in the firing rate variability of pairs of cortical neurons (Gawne and Richmond, 1993; Zohary et al., 1994; Lee et al., 1998). The contribution of correlated discharge has largely been dismissed as reflecting coding redundancy in a noisy system. However, correlations in response variability can provide additional information that is available only from the assembly (Deadwyler and Hampson, 1996;Hatsopoulos et al., 1998b; Oram et al., 1998). In the present study we demonstrate that statistical dependencies in the trial-to-trial firing rates of pairs of neurons can be an additional source of information about motor behavior.

## MATERIALS AND METHODS

*Task.* Two male monkeys *(macaca fasicularis)* were operantly conditioned to perform an instructed-delay task consisting of visually guided planar reaching movements from a central holding position to a radially located target position (Georgopoulos et al., 1986; Hatsopoulos et al., 1998a). Animals moved a two-joint manipulandum in the horizontal plane to direct a cursor from a central hold position to one of eight possible radially positioned targets (1 cm^{2}, 6 cm away from start position) viewed on a computer monitor. A trial was composed of three epochs: a hold period during which time the monkey had to maintain the cursor at the hold position for 0.5 sec, a random 1–1.5 sec “instructed delay” period during which the target for the forthcoming action appeared but movement was withheld, and a “go” period initiated by target blinking (mean reaction time, ∼365 msec). Manipulandum position was monitored using a digitizing tablet (Numonics) sampled at 72 Hz. Monkeys generally completed >30 movements in each of the eight directions during a recording session.

*Array recordings.* When >85% correct trials were attained, an array of 100 microelectrodes was implanted in the MI arm representation, medial to the spur of the arcuate sulcus in the mediolateral dimension, abutting the central sulcus in the anteropoterior dimension. The microelectrode arrays (Bionic Technologies, Salt Lake City, UT) consisted of 100, 1.5-mm-long platinized tip, silicon probes arranged in a square grid (400 μm on center) with probe impedances between 200 and 500 kΩ (1 nA, 1 kHz sine wave; Nordhausen et al., 1996). Twenty-two in the central matrix of the 100 available electrodes were selected in advance for electrical connection to one of two percutaneous connectors using 25.4 μm (1 mil) Teflon-coated platinum–iridium (10%) wire (10IR1T; MedWire, New York, NY). An extra wire (50.8 μm diameter, Pt-Ir 20%) with ∼5.0 mm of the terminal insulation removed was inserted subdurally and used as a reference. For added ease of handling, the 12 wires were held together in four or five spots with silicone elastomer (MDX4–4210; Dow Corning, Midland, MI). Each percutaneous connector consisted of a 12-pin connector (FR-12S-6; Microtech, Pittsburgh, PA) set in a custom-designed titanium pedestal. The back of the array and percutaneous connectors were coated with silicone elastomer to mechanically protect the wires and maintain electrical insulation at the bond pad sites.

Animal care and surgical procedures used in the identification of the arm area in monkey MI are described elsewhere (Donoghue et al., 1998), as are the procedures for implanting the array and closing the surgical defect (Rousche and Normann, 1992). All surgical procedures conducted for this study were in accordance with Brown University Institutional Animal Care and Use Committee-approved protocols and the Guide for the Care and Use of Laboratory Animals (National Institute of Health publication no. 85-23, revised 1985).

Recordings began after a 1 week period of postsurgical recovery. Custom shielded cabling carried signals from either 11 or 22 electrodes to a custom 16-channel amplifier (gain, 25,000; −3 bB points 250 Hz and 7.5 kHz; 5th order Butterworth anti-aliasing filter; input impedance, 10^{9} Ω). Signals were simultaneously digitized at 20 kHz (Experimenter’s Workbench; DataWave Technologies, Longmont, CO) and processed to extract action potentials. All candidate action potential waveforms (−0.5 to +1.5 msec around threshold crossing) exceeding ∼2× background noise were stored for off-line processing and analysis. A given data set included the simultaneously recorded neural responses from 16 electrodes and kinematic data in the form of the digitized x-y position of the manipulandum.

Array electrodes often recorded action potentials originating from more than one neuron (Nordhausen et al., 1996; Maynard et al., 1997) and ranged from multiple nonclassifiable units to well isolated potentials. Classification of the action potential waveforms was accomplished using a semiautomatic spike-sorting algorithm based on eight waveform parameters (Autocut; DataWave Technologies). The classification templates resulting from this process were subsequently verified and adjusted manually. In a small number of the identified neurons, it was likely that a few appropriate waveforms were excluded in the process of removing incorrectly classified waveforms. Units with stable, identifiable waveform shapes and a refractory period >1 msec demonstrated in autocorrelograms were considered to be well-isolated single units. Action potentials were reduced to an event time series for further analysis. Based on the lack of a clear refractory period, 13 of the 85 (15.3%) recorded neurons were excluded from further analysis. Some of the remaining neurons showing a few events in the refractory period were nevertheless included in the analyses to maximize sample size and labeled as multiunit neurons (Table1). As will be shown below, results from data sets that included these multiunits were similar to those sets with only single units.

For each identified unit, the mean firing rate during each movement trial was calculated from the number of action potentials present in the 600 msec interval beginning 200 msec before movement onset and ending 400 msec afterwards. This interval most commonly showed the greatest firing rate change from background for the recorded units. Movement onset was defined as the time that hand tangential velocity exceeded 1.0 cm/sec.

*Statistical analyses.* The statistical significance of an observed correlation coefficient for a pair of neurons was assessed by comparing it with a correlation coefficient computed after the sequence of trials associated with one neuron was randomly shuffled relative to that of the other neuron. Shuffling removes trial specific correlations and creates an “independent” condition similar to that which would be obtained in serially recorded (i.e., one neuron at a time) data. For each cell pair in each movement direction, 10 unshuffled subsets of data (>29 movement trials per subset) were selected from the complete data set and used to compute 10 correlation coefficients and then averaged. The temporal sequence of these trials was then randomly shuffled, and a set of 10 control correlation coefficients was computed. The observed average correlation coefficient for each pair of cells in each of the eight movement directions was assessed for significance by comparing it to the mean of the 10 control correlation coefficients using a two-sample *t* test (*p* < 0.01).

In practice, the covariance matrices that characterized the second-order multivariate Gaussian model used for classification (see Fig. 8) for more than two neurons were not always invertible because singularities result from small numbers of trials. To avoid this problem, we resorted to generating realizations of a multivariate Gaussian distribution with a covariance structure (and mean) given by the actual neural data using a pseudorandom number generator without having to invert the covariance matrix (Johnson, 1994). One hundred realizations were generated for each of the eight movement directions. Each trial of the actual neural data were then compared to the*k*-nearest neighbor realizations (using a Euclidean distance measure) and assigned to the movement direction that generated the majority of the *k*-nearest neighbor realizations. The nearest neighbor classifier error rate never exceeds twice that of the optimal Bayesian classifier and will approach the Bayes’ error rate as the number of realizations and *k* go to infinity (Duda and Hart, 1973). A variety of values for *k* were used and were found to generate qualitatively similar results. The results reported here were based on a value of 10. To ensure that the classifier was not overfitting the neural data used to construct it, the classifier was cross-validated by testing it on neural data that was not used to build it. Cross-validation was implemented by leaving one trial out and generating the mean and covariance parameters based on the remaining trials. This was repeated for all trials in the data set.

## RESULTS

The data used for this study were obtained from five recording sessions in two monkeys made over a 1 year period. Most neurons demonstrated a marked modulation in their average firing rate around arm movement but did not modulate their firing rate strongly if at all during the instructed-delay period (Fig.1*A*). To characterize the directional tuning of each neuron, a linear regression was used to fit a cosine function to the observed firing rates in the eight movement directions (Schwartz et al., 1988). The directional tuning curves of 64 of the 72 (89%) units were fit by a cosine function (*p* < 0.05). The peak of this cosine function identified the preferred direction of that cell (Schwartz et al., 1988). The aggregate of preferred directions did not uniformly represent the movement space (Fig. 1*B*) in either monkey, particularly between 0 and 90°, over all five data sets (Fig.1*B*, *bottom plot*). This nonuniformity is most likely attributable to limited sample size because others have found directional tuning spanning the entire two-dimensional workspace (Georgopoulos et al., 1986).

We recorded data from the same location in one animal in sessions separated by 23, 26, and 329 days. A preliminary stability analysis based on the inspection of waveforms, autocorrelograms, and perievent histograms in one monkey suggested that >50% of the sites recorded action potentials from different neurons after a 1 month period and 100% after 14 months. Therefore, some of the recording sessions may have contained partially overlapping populations.

### Variability in neural discharge

As is evident in Figure 1*A*, the modulation in MI firing during repeated movements in the same direction vary somewhat from trial to trial. Response variability for the set of MI neurons was quantified by comparing the SD of the observed firing rates for all of the movements performed in a given direction to the mean firing rate for that direction. The relationship between the trial-to-trial variability of the response and the mean firing rate is best fit by a power function (Fig. 2). The dependence of the variability (SD) on the mean firing rate (M) is described by a power function of the form: SD = 1.44 M^{0.50}(*r*^{2} = 0.85). This relationship is consistent with the findings of Lee et al. (1998) who derived a relationship of SD = 1.21 M^{0.59}(*r*^{2} = 0.77). The variance exceeds that predicted by a Poisson count in which the mean equals the variance. Although the origin of this excess variance is unclear, we will show that some of this variability as revealed from pairwise correlations can be accounted for by the covariation between cells.

### Covarying neural discharge

To study correlations among MI neurons, the following equation (Eq. 1) was used to normalize responses to their respective*z*-scores on a trial-by-trial basis (Crow et al., 1960; Lee et al., 1998):
Equation 1In Equation 1, *i* indexes a particular neuron,*k* specifies the trial number, *d* specifies the movement direction, *m*^{d} is the mean, and ς^{d} is the SD of the firing rates for movements in a given direction, *d*.

Although responses of single neurons were variable from trial to trial, this variability was correlated for many pairs of MI neurons. Figure6*A* (*left*) depicts an exceptional example of correlated activity on a trial-to-trial basis for a pair of neurons simultaneously recorded from MI. The correlation strength across MI neuron pairs was evaluated by computing the correlation coefficient,*ρ*_{ij}^{d}*,*between the responses of neuron *i* and *j* for each movement direction *d*, separately:
Equation 2where *E*[… ] represents the expected value of the product of the firing rate *z* scores for neurons *i*and *j* for movement direction *d*. This expected value was estimated from the sample mean over all trials for one movement direction. A Fisher’s *z* transform was used to convert the correlation coefficients so that they became normally distributed and so that populations of correlation coefficients could be more easily compared (Lee et al., 1998):
Equation 3Correlation coefficients were computed for each pair of cells, separately for each direction. The average correlation for all pairs of cells in the eight movement directions was 0.21 ± 0.34. The slight bias toward positive correlation in the trial-to-trial variability is consistent with other reported correlations in cortical neurons [e.g., 0.12 by Zohary et al. (1994) and Lee et al. (1998)]. In the entire population of simultaneously recorded MI neurons for each of the eight movement directions, significant correlations (*t* test; *p* < 0.01) were found in 2538 of the 3976 (63.8%) pairs. Figure 3 shows the broad distribution of the 2502 significant correlation coefficients between pairs of simultaneously recorded MI neurons. Note that a few of the interactions present in the pairs of MI neurons for movements in a particular direction were quite strong (i.e., correlation coefficients of ∼0.8 or higher).

### Correlation and directional preference

The dependence of the pairwise response correlation on direction preference was investigated. A strong relationship between the similarity of directional tuning and response correlation would be expected if the correlation resulted from redundant information about movement for cell pairs. Directional tuning similarity for each pair of simultaneously recorded MI neurons was evaluated by computing the dot product between the preferred directions of the two cells. This procedure gives a quantitative measure of similarity between two vector quantities such as direction because the dot product is proportional to the cosine of the angle subtended by the two vectors. Thus, pairs of cells with a dot product near 1.0 have very similar preferred directions. Correlation coefficients between the responses of neuron pairs as a function of the similarity in their directional tuning curves are shown in Figure 4. Although a linear (slope = 0.09; *r*^{2} = 0.74) relationship was found between the mean correlation coefficient and the similarity of the preferred directions, the variability in correlation strength was very large over different cell pairs. That is, different cell pairs could exhibit very different correlation coefficients even when the similarity in their preferred directions was common. Thus, there appears to be a weak relationship between the directional preference and correlation strength when correlations are measured on a broad time scale (i.e., 600 msec).

### Correlation and interelectrode distance

Based on the distribution of MI afferents from other regions of cortex, cells that are close to each other have a higher likelihood of receiving common inputs, thus modulation of the shared input could result in the observed correlation in firing of neuronal pairs. Shared input has been a common explanation for coding redundancy when cell pairs are recorded on single electrodes (Gawne and Richmond, 1993). The possibility that correlations between MI neurons were the result of similar covarying inputs was investigated by plotting the magnitude of the correlation coefficients as a function of the cortical distance between the electrodes on which the neurons were recorded (Fig.5). For the significant correlations in our population of simultaneously recorded neurons, there was no significant relationship between the strength of the correlations and the cortical distance between the neurons (linear regression,*p* < 0.36). As is evident in the plot, pairs of neurons recorded from the same electrode (distance = 0) were no more likely to be correlated than neurons separated by distances as much as 2 mm. This result suggests that these correlations emerge from spatially broad interactions and not from local common input.

### Directional variations in correlation

Although correlated activity between cortical neurons has been previously observed, we discovered that correlation coefficients between pairs of MI neurons could vary with movement direction [seeHatsopoulos et al. (1998b), for directional changes in synchrony measured on a fine time scale]. Notice how the same cell pair exhibits strong correlated activity on a trial-by-trial basis for one movement direction (Fig. 6*A, top*) but much weaker correlated activity for a different movement direction (Fig. 6*A, bottom*). Figure6*B* provides an example of directional variations in correlation between two MI neurons over all eight directions. As is evident from the firing rate tuning curves of the two cells, they have similar preferred directions of ∼135°; yet they engage in significant correlated firing for movements in the 45 and 270° directions. A one-way ANOVA for direction revealed that 78% (387/497) of the cell pairs exhibited significant differences in correlation with movement direction (*p* < 0.05). This finding indicates that correlation strength is not necessarily a fixed property of neuronal pairs but rather can depend on movement direction.

### Covarying neural discharge and directional coding

We evaluated whether correlations could contribute to directional coding by comparing the performance of statistical models that incorporate second-order relations with those that assume independence of neurons in the population. These statistical models were used to classify single trials of neuronal data into one of eight classes corresponding to each of the eight movement directions. Bayesian estimation techniques (Sanger, 1996; Oram et al., 1998; Zhang et al., 1998) as well as approximations to Bayesian estimation were used to determine the direction of the arm movement. Bayesian techniques use probabilistic descriptions of the relationship between the firing rate of a neuron and the direction of arm movement to return the probability of a movement direction given the firing characteristics of a population of cells within a trial. This probabilistic description was constructed by assuming that the number of action potentials across all the neurons in the population,
for movements in a particular direction, *d,* followed a multivariate Gaussian distribution:
Equation 4where *N* is the number of neurons in the population,
is the mean population response, and *C*^{d} is the covariance matrix for movements in direction, *d*. The vertical lines surrounding the covariance matrix denote its determinant. The superscript *T* denotes the transpose operator. The population response on a given trial,
, and the mean population response,
^{d}, were obtained directly from the experimental data.

To determine the direction of the arm movement from the response of the ensemble of MI neurons, Bayes’ rule was used to transform the probability of observing a response in a particular movement direction,*P*(
‖*d*), into a probability that a movement was performed in a particular direction given a particular pattern of neural firing (i.e., the neural “response”):
Equation 5where *P(d)* is the prior probability of a particular movement direction and*P*(
) is the prior probability of a population firing pattern.

For each trial, these eight probabilities were collapsed into a single estimate for the direction of movement by choosing the direction with the highest probability of occurrence (i.e., the maximum *a posteriori* estimate or MAP estimator). If the direction chosen was the same as the actual direction in which the monkey moved, then the estimation was successful (chance, 12.5%). Because of the limited data available for these analyses, cross-validation of the model was performed by removing each trial from the data set, forming the various probability distributions on the remaining data, and estimating the direction of movement for the firing rates in the removed trial.

The role that the response correlations might play in directional coding was tested by comparing two versions of the Gaussian model defined in Equation 4. The first model assumed that the responses of each of the neurons were conditionally independent. This was accomplished by setting the off-diagonal terms of the covariance matrix, *C*^{d}, to zero and only estimating on-diagonal terms, i.e., the variances of each of the neurons. In other words, the responses of all cells were assumed to be uncorrelated with each other and, therefore, independent of each other given that they are normally distributed. It is important to stress that this is a conditional independence assumption. That is, conditional on a particular movement direction, the probability that one neuron gives a particular response is independent of the probability that any other neuron gives a particular response. Clearly, cells with similar preferred directions are unconditionally dependent. The second model, on the other hand, incorporated the observed second-order correlations by estimating the off-diagonal elements of the covariance matrix. To illustrate the difference between the two models, data based on the activity of two neurons were classified into one of two possible movement directions using either the independent model (Fig. 7*A*) or the second-order model (Fig. 7*B*). Trials from one data set are correctly classified (solid symbols) or incorrectly classified (hollow symbols) as 180° (circles) or 270° (triangles) movements. In addition, contour plots show the isoprobability lines representing the joint probability of activity of two cells conditional on either the 180° direction (dashed lines) or 270° direction (sold lines). Notice how the second-order model (Fig. 7*B*) fits the data better resulting in fewer incorrectly classified trials.

The classification results for the independent and second-order models based on the complete ensembles of recorded cells (instead of just pairs of neurons) are shown in Figure8*A*. The shaded bars in the figure show the success of the model that treated the neurons as independent encoders of movement direction in determining the direction of movement from a population of simultaneously recorded MI neurons. The dark bars show the success of the second-order model that incorporated the pairwise covariances in predicting the direction of movement from a simultaneously recorded population of MI neurons. In all cases, incorporating the second-order interactions between neurons improved classification of movement direction (Fig.8*A*). On average, incorporating the covariance structure resulted in correct determination of movement direction in 72% ± 17% (±1 SD) of the trials; the model without the correlations correctly estimated the direction of movement in 61% ± 12% of the trials. The difference between these two populations (11%) was significant (paired *t* test, *t*(4) = 3.97, one-tailed, *p* < 0.05); this test demonstrates that these correlations provide additional information about direction not present in the firing rates of individual neurons treated as independent coding elements. In both cases, the models performed significantly better than the chance rate of 12.5%.

We next investigated to what extent the observed changes in correlation with direction could account for the 11% improvement in the performance of the second-order model. If we treat the correlation strength as composed of a fixed DC component and a variable AC component, the bulk of the improvement in directional estimation may be accounted for by incorporating the fixed correlational structure into the second-order model. On the other hand, the observed changes in correlation with movement direction may be contributing substantially to the improvement of the second-order model. To tease apart these two contributions, we compared the performance of the original second-order model having variable correlations to a second-order model whose covariance matrices were fixed to the average of their values over all eight movement directions. Figure 8*B* shows that the original second-order model does consistently better than the model with a fixed covariance structure. The model with fixed covariance performed only 5% better than the independent model compared to an 11% improvement with the original second-order model. Thus, 55% of the improvement [(11 − 5)/11] in performance of the second-order model is because of variations in correlation with movement direction.

### The nature of information contained in covarying neural discharge

Correlations between the responses of neurons can either encode information redundantly, independently, or synergistically (Gawne and Richmond, 1993). When information is redundantly encoded in a pair of neurons, there is less information in the simultaneous firing rates of the two cells than if the two neurons independently represented information (i.e., if one neuron indicates a leftward direction with 100% reliability, having a second neuron that indicates left is redundant). Likewise, when two neurons synergistically encode information, the simultaneously recorded firing rates of the two neurons contain more information than if the neurons were taken independently (i.e., the coactivation of a left direction and a right direction neuron may indicate a new movement direction). Because accurate estimates of information content in rate codes requires more data than we were able to collect (Golomb et al., 1997), we used statistical models that were constrained by the experimental data to determine the nature of information available from pairs of neurons.

The model was constructed using (1) the cosine directional tuning functions for neurons whose firing rate modulation was well-described by a cosine, (2) a model for the dependence of the response variability of a single cell on its mean firing rate (Fig. 2), and (3) the observed correlations in the responses of simultaneously recorded cell pairs. These three parameters were derived from the data and then used to create discrete bivariate Gaussian models,*P(r*‖*d),* for pairs of neurons as well univariate Gaussian models for individual neurons. These models were similar to those used in the Bayesian classification analysis described above except that they used discrete instead of continuous normal distributions. Using Equations 5 and 6, it was possible to compute the average mutual information *I*(*D*;*R*) between the joint response of the pair of neurons and the direction of the arm movement. It was compared to the sum of the information from the two constituent neurons (Gawne and Richmond, 1993). If the correlations represent directional information redundantly, the ratio of the joint mutual information to the sum of the mutual information from the two constituent neurons treated separately is <1.0. If the ratio is >1.0, the correlations between the simultaneously recorded neurons contribute more information than is available from the two cells treated independently.
Equation 6where the angled brackets represent an average over neuronal responses. Figure 9 shows the distribution of the information fraction for 342 pairs of neurons. This distribution indicates that covarying discharges in a pair of MI neurons can synergistically encode information about movement direction, although for many cases the coding is redundant. The significant shift in the mean of this distribution (1.064 ± 0.006; *p* < 0.001; Wilcoxon signed-rank test) suggests that that correlations in the trial-to-trial response variability of MI pairs as a group are more likely to increase than decrease information content. In some cases, the presence of correlations can increase the representational capacity of the pair of neurons in excess of 30% over the independent case.

## DISCUSSION

The present study demonstrates the capability of neuronal populations to represent information in their joint activity. The second-order covariations in firing between MI neurons provide directional information above that found in single neuron firing rates. Directional classification of data on a trial-to-trial basis was consistently better when these interactions were considered. Information-theoretic analysis demonstrated that the joint activity of certain neuron pairs could provide >30% additional information about direction (Fig. 9). These findings suggest that higher-order ensemble codes that incorporate the interactions of a population of neurons can make important contributions to cortical representations [seeHatsopoulos et al. (1998a) for a comparison between ensemble and standard population codes] if these codes can be exploited. We have not demonstrated that these codes are actually used but that additional information is available when these interactions are considered.

### Noise in single neuron firing rates

Although it has been shown that single MI neurons carry information about movement direction in their average firing rate (Georgopoulos et al., 1982), the firing rate variability is higher than would be expected from a random Poisson process [Fig. 2; see also Fig. 2 in Lee et al. (1998)]. However, we know that precisely timed and controlled movements can be generated despite this variability. To overcome the problem of random noise, precise movements could be represented by averaging the responses of a population of MI neurons (Georgopoulos et al., 1986). When noise sources are statistically independent across different neurons, the effect of noise on the population average decreases by a factor of one over the square root of the number of neurons, thereby increasing the signal-to-noise ratio. This explanation becomes less tenable as noise becomes more correlated across cells (Zohary et al., 1994). We have shown that this noise is partially correlated, in agreement with others (van Kan et al., 1985;Gawne and Richmond, 1993; Zohary et al., 1994; Lee et al., 1998), but unlike these studies we have emphasized that this correlated noise is not fixed but rather varies with behavioral condition, i.e., movement direction.

Our approach considered correlated noise as part of the signal and demonstrated that it carried additional information about movement direction above that obtained from the firing rates of individual neurons acting as independent encoders. Thus, rather than limiting the information capabilities of neural populations, correlations may be able to extend them as an additional type of signal. A key feature of this approach is the finding that the firing covariation between neurons varies across different movement directions. As Figure8*B* demonstrates, a major contribution to the improvement in classification performance of the second-order model was because of the fact that the directional variations in firing covariance among motor cortical neurons were incorporated into the model. That is, the firing covariances were preserved for each movement direction unlike other studies, which pooled the covariances over different stimulus or movement directions (Zohary et al., 1994; Lee et al., 1998).

These findings complement other results from MI that demonstrate that information can be obtained from correlated activity between neurons on a finer time scale. Riehle et al. (1997) have shown that synchronous pairwise firing on a millisecond time scale can signal the expectation of visual cues instructing the monkey to initiate a movement in the absence of rate modulation. In addition, results from our laboratory have demonstrated that directional information is available in fine synchronous interactions of neurons that cannot be accounted for by independent rate modulation (Hatsopoulos et al., 1998b). We have not fully examined the relationship, if any, between broad correlations in spike count measured on the time scale of hundreds of milliseconds and fine synchrony measured on the millisecond time scale. However, preliminary results suggest that they are relatively independent phenomena. The flat relationship between broad correlation and interelectrode distance (Fig. 5) is inconsistent with the relatively strong dependence of fine synchrony with distance that has been observed previously (Kwan et al., 1987; Hatsopoulos et al., 1998b, their Fig. 1*D*). This also suggests that the two processes may be independent.

Our results suggest that the observed firing variability may be considered a potential signal when viewed in the context of an ensemble of neurons. This may be true not only in MI but throughout cortex. Studies in the hippocampus and somatosensory cortex using different ensemble-based approaches have also found that classification of behavior or sensory localization improves when higher dimensional features of a neural ensemble are incorporated (Deadwyler et al., 1996;Nicolelis et al., 1998).

### Source of correlation

Although correlated activity between cortical neurons has been widely detected, the source of these interactions has not been established. One possibility is that correlations are a consequence of common input from other neurons. Nearby neurons, which have been the focus of most recording studies, might be expected to be correlated because they would share thalamic projections (in MI, Shinoda and Kakei, 1989). However, the bulk of our neurons would not likely be driven by common thalamic sources because our electrodes were spaced by at least 400 μm, which is a considerable distance for two cells to receive strong, common input. Moreover, the observation that the correlation strength does not vary with distance argues further that common input is not a major source of correlation.

An alternative hypothesis is that correlations arise from large-scale reciprocal connections among cortical neurons (Juergens and Eckhorn,1997). Neurons in motor cortex form a network of connections with each other with a high degree of lateral interconnectivity (Ichikawa et al., 1985; Huntley and Jones, 1991; Hess and Donoghue, 1994). Moreover, these lateral connections can become stronger with the acquisition of new motor behaviors (Rioult-Pedotti et al., 1998). Functional studies have also shown that a large, distributed area of the MI arm area activates when simple reaching movements are made, suggesting that motor cortex is not organized as a set of independent columns (Sanes et al., 1995).

### Functional role of correlations

It is unclear what functional role these correlations may be playing, if any. One possibility is that they are the consequence of anatomical connections between cortical neurons that may be important in promoting population codes. Recent theoretical work has suggested that lateral connections could be useful in refining a coarse code among a population of cortical neurons (Pouget et al., 1998). A particular arrangement of lateral connections between directionally selective neurons would reduce the effects of neural noise by creating a smooth “hill” of activity over the set of cortical neurons. According to the neural model they develop, the correlation between two directionally selective neurons is predicted to be proportional to the product of the derivative of their tuning curves. Although this prediction has not been systematically tested with our data, the variations in correlation demonstrated in Figure 6*B*are somewhat consistent with it.

Another hypothesis specific to motor cortex is that these correlations are used to create temporary synergies among muscle groups. Layer 5 pyramidal neurons in primary motor cortex form part of the corticospinal tract and, thus, directly and indirectly synapse onto motor neuron pools that innervate skeletal muscles. Hepp-Raymond et al. (1996) have recently made electromyographic recordings from hand muscles that are active during precision grip movements. It was observed that particular muscles became correlated and, more importantly, the correlation strength varied with the overall level of force. These variable muscle synergies may be created by linking different neurons in motor cortex.

### Methodological limitations

This work has demonstrated that the presence of second-order interactions in a population of MI neurons can significantly add to the representation of directional information. This is an initial step toward examining higher-order relationships among ensembles of neurons. We restricted ourselves to a relatively simple multivariate model (i.e., the Gaussian model) that considered only pairwise correlations between neurons. This was, in part, because of the limited amount of data that could be acquired from a behaving animal and is a major challenge in the endeavor of discovering neural ensemble representations. It is unknown, however, whether even higher-order moments may account for more of the response variability and, therefore, provide more information about direction or any other parameter. As multielectrode recording becomes more prevalent, novel statistical methods will be needed that are mathematically sound and tractable, and which incorporate multidimensional interactions.

Although we have shown that higher-order representations provide more information about movement direction, these experiments do not demonstrate that these higher-order representations are actually being used by the animal to control motor behavior. In fact, most if not all studies that postulate a coding scheme based on the firing rates of single neurons have the same problem. A recent study of olfactory coding in an invertebrate preparation has attempted to show that correlated neural activity as revealed by field potential oscillations is necessary for fine olfactory discrimination (Stopfer et al., 1997). This was accomplished by pharmacologically desynchronizing a population of neurons without affecting their individual firing properties and showing that olfactory discrimination was impaired. In the future, similar perturbation experiments could be done in the cortex of mammals to determine whether higher-order interactions among neurons are necessary to guide behavior.

## Appendix

The probability of observing a number of action potentials across an ensemble of *N* motor cortical neurons,
, given the monkey’s movement direction, *d*, was modeled as a multivariate Gaussian distribution:
Equation 4Eight of these distributions were specified, one for each direction. The parameters that fully specify each distribution are: (1) the mean spike count over the ensemble of neurons for a particular movement direction,
^{d}. This vector of mean values is *N*-dimensional; (2) the covariance matrix, *C*^{d}, computed for a particular movement direction. This matrix is *N × N*. An element of the covariance matrix is computed as follows:
where *M*^{d} is the number of trials for direction *d*. The inverse of the covariance matrix is defined as follows:*C*^{d}*C*^{d−1}where *I* is the identity matrix. The vertical lines in the denominator of the above equation represents the determinant of the matrix.

The mutual information, *I(D*;*R)*, between the movement direction, *D*, and the response of the neural ensemble, *R*, is defined as follows:
Equation 6The *a priori* probability of a movement direction,*P(d)*, is determined by the experimenter and was specified to be uniform. The *a posteriori* probability of a movement direction given the response of the ensemble,*P*(*d*‖
), was estimated by using Bayes’ rule (Eq. 5). This equation can be decomposed into three parts: (1) the logarithm of the ratio of the*a posteriori* to *a priori* probability distributions (the term in the brackets). If the firing rates of the ensemble do not alter one’s estimation of the probability of a particular direction, then the ensemble response doesn’t provide information about direction. In this case, the ratio of distributions in Equation 6 will be 1, and the log of 1 is zero, which confirms our intuition. On the other hand, if the ensemble response does alter one’s estimation of the probability of a particular direction, then the ratio is not 1. If the *a posteriori* probability is large, the ratio will be >1, which will give a positive log. (2) The expectation or average of the information over all movement directions conditional on the ensemble response (the sum sign and its terms). The information is weighted more for directions that are conditionally more likely. (3) The average over all ensemble responses (is represented with the angular brackets).

## Footnotes

This work was supported by National Institutes of Health Grants NS25074 and NS10269-02, National Science Foundation Grant SBR-N720368, the Utah Center of Excellence in Neural Interfaces, and the W. M. Keck Foundation. E.M.M., N.G.H., and J.P.D. are responsible for the concept, analysis, and writing of the this manuscript. All other authors were involved in the surgical and technical aspects of this work, as well as in the data acquisition. We thank Elie Bienenstock and Stu Geman for helpful comments on this manuscript. We also thank Beth Travers and Gerhard Friehs for their technical assistance.

The first two authors are the primary authors of this paper.

Correspondence should be addressed to Nicholas Hatsopoulos, Box 1953, Department of Neuroscience, Brown University, Providence, RI 02912.