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 armreaching task. Pairs of simultaneously recorded cells revealed significant correlations in their trialtotrial 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. Informationtheoretic 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; RioultPedotti 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 armreaching 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 trialtotrial 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 instructeddelay 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 twojoint 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.5mmlong 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). Twentytwo 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) Tefloncoated platinum–iridium (10%) wire (10IR1T; MedWire, New York, NY). An extra wire (50.8 μm diameter, PtIr 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 12pin connector (FR12S6; Microtech, Pittsburgh, PA) set in a customdesigned 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 Committeeapproved protocols and the Guide for the Care and Use of Laboratory Animals (National Institute of Health publication no. 8523, 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 16channel amplifier (gain, 25,000; −3 bB points 250 Hz and 7.5 kHz; 5th order Butterworth antialiasing 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 offline 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 xy 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 spikesorting 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 wellisolated 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 twosample t test (p < 0.01).
In practice, the covariance matrices that characterized the secondorder 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 theknearest neighbor realizations (using a Euclidean distance measure) and assigned to the movement direction that generated the majority of the knearest 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 crossvalidated by testing it on neural data that was not used to build it. Crossvalidation 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 instructeddelay 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 twodimensional 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 trialtotrial 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 respectivezscores on a trialbytrial basis (Crow et al., 1960; Lee et al., 1998):
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 trialtotrial 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:
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 trialbytrial 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 oneway 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 secondorder 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,
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(
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, crossvalidation 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 offdiagonal terms of the covariance matrix, C ^{d}, to zero and only estimating ondiagonal 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 secondorder correlations by estimating the offdiagonal 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 secondorder 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 secondorder model (Fig. 7 B) fits the data better resulting in fewer incorrectly classified trials.
The classification results for the independent and secondorder 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 secondorder 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 secondorder 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, onetailed, 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 secondorder 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 secondorder model. On the other hand, the observed changes in correlation with movement direction may be contributing substantially to the improvement of the secondorder model. To tease apart these two contributions, we compared the performance of the original secondorder model having variable correlations to a secondorder model whose covariance matrices were fixed to the average of their values over all eight movement directions. Figure 8 B shows that the original secondorder 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 secondorder model. Thus, 55% of the improvement [(11 − 5)/11] in performance of the secondorder 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 welldescribed 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.
DISCUSSION
The present study demonstrates the capability of neuronal populations to represent information in their joint activity. The secondorder covariations in firing between MI neurons provide directional information above that found in single neuron firing rates. Directional classification of data on a trialtotrial basis was consistently better when these interactions were considered. Informationtheoretic analysis demonstrated that the joint activity of certain neuron pairs could provide >30% additional information about direction (Fig. 9). These findings suggest that higherorder 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 signaltonoise 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 secondorder 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. 1D). 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 ensemblebased 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 largescale 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 (RioultPedotti 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 Bare 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. HeppRaymond 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 secondorder interactions in a population of MI neurons can significantly add to the representation of directional information. This is an initial step toward examining higherorder 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 higherorder 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 higherorder representations provide more information about movement direction, these experiments do not demonstrate that these higherorder 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 higherorder 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,
The mutual information, I(D;R), between the movement direction, D, and the response of the neural ensemble, R, is defined as follows:
Footnotes

This work was supported by National Institutes of Health Grants NS25074 and NS1026902, National Science Foundation Grant SBRN720368, 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.