we analyzed 12 neuronal datasets recorded from two human clinical trial participants with tetraplegia ( hs1 and hs3 ) and four monkeys ( mla , mcl , mco and mab ) . these included m1 recordings taken while human participants performed two sessions of a  neural cursor center - out task ( that is , a task where the participant used , via a neural interface , his or her own recorded m1 spiking activity to move a computer cursor to targets radially positioned on the computer screen ) , m1 recordings from two monkeys ( mla and mcl ) , each performing two sessions of a task requiring planar point - to - point reaches toward targets randomly placed in the workspace , simultaneous m1 and pmv recordings from a monkey ( mco ) performing two sessions of a task that required reach and grasp toward moving objects in a three - dimensional workspace , and simultaneous m1 and 5d recordings from a monkey ( mab ) performing two sessions of a pursuit tracking task that required planar hand movements . minimum inter - electrode distance in the 10  10 microelectrode array was 400 m . maximum inter - electrode distance , including electrodes in two arrays , was ~2 cm ( supplementary fig . 1 ) . to assess the predictive power of spiking histories , we first computed the probability that any given i neuron xi , t was going to spike at time t , conditioned on spiking histories t from the past 100 ms up to ( but not including ) time t. without further constraints , direct estimation of conditional probability distributions for high - dimensional systems is typically an intractable problem , leading to combinatorial explosion and requiring amounts of data that grow exponentially with the number of neurons in the ensemble . instead , we took advantage of the fact that this conditional probability can be computed as16 
 equation ( 1)pr(xi , t=1|t)t(t|t) where i(t| t ) is the instantaneous spiking rate ( conditional intensity function ) of the neuron and  = 1 ms in our discrete time representation , and used a simplified model to capture the relationship between the instantaneous rate and spiking histories equation ( 2)log(i(t|t))=i+k1,ixi+jik2,i , jxj 
  the term i relates to a background level of spiking activity , xi is the spiking history in the specified time interval for the i neuron , i = 1 , 2 ,  , n recorded neurons , and k1,i and k2,i , j denote temporal filters related to intrinsic and ensemble history effects , respectively . once an instantaneous spiking rate model was fitted , the estimated probability of a spike at any given time bin , conditioned on intrinsic and ensemble spiking histories , was easily computed using equation ( 1 ) ( the spike prediction approach is shown using cell 34a ( hs3 ) as the example target neuron ; fig . 1 ) . on the basis of the estimated conditional spiking probabilities , we used a standard tool , the roc curve analysis19 , to assess the predictive power of spiking history models ( see online methods ) . spike prediction in 1-ms time bins based on spiking histories was substantial . intra - areal ensemble history  based predictions in m1 resulted in high true - positive rates while maintaining low numbers of false positives . for example , it was possible to correctly predict 80% of spikes in neuron 34a ( participant hs3 ) with a false - positive rate of less than 5% ( fig . similarly , an 80% true - positive rate was obtained for neuron 16a ( monkey mla ) with a less than 10% false - positive rate ( fig . this predictability was substantial despite the weak pair - wise correlations that we observed among all neuronal pairs in the recorded datasets for hs3 and mla ( fig . 2c , we obtained a more comprehensive assessment of predictive performance by computing the area under the roc curve . the area under the curve ( auc ) is a global summary statistic ; that is , it depends on both the true- and false - positive rates and on all of the possible thresholds on the spiking probability . the auc gives the probability that , when two samples are randomly drawn from the data ( one containing a spike , the other not ) , the conditional intensity model will assign a higher probability ( i.e. a higher instantaneous spiking rate ) to the sample with a spike20 . it therefore provides an assessment of the discriminatory power for predictive variables under a given model . it approaches 0.5 for a chance level predictor , that is , a predictor having false- and true - positive rates along the diagonal , and it equals 1 for a perfect predictor . 2a ) . we further corrected the auc by subtracting chance - level predictions estimated from the actual data and scaled it to obtain a quantity that ranged from 0 ( no predictive power ) to 1 ( perfect prediction ) . we refer to this corrected and normalized auc measure as the predictive power ( see online methods ) . more than 50% of examined neurons in monkey mla showed a predictive power higher than 0.5 ( median = 0.55 , range = 0.10.89 ) ( fig . 3 ) . prediction in subjects mcl and hs3 showed similar performances , with medians of 0.45 and 0.49 and ranges of 0.180.91 and 0.230.97 , respectively . predictive power was lower for participant hs1 ( median = 0.30 , range = 0140.69 ) . we conjectured that this was probably a result of the smaller number of neurons ( n = 21 , n = 22 ) that were recorded and used for prediction in this participant . to test this conjecture , we computed the predictive power on the basis of several smaller ensemble subsets ( n = 22 neurons ) randomly chosen out of the 110 neurons from hs3 s dataset ( session 2 ) . data from hs3 provided a good reference , as neurons were recorded from m1 under the same task condition as for hs1 . the distribution of predictive power values , based on this collection of smaller ensemble subsets , was similar to the one obtained for hs1 , supporting our hypothesis ( supplementary fig . 2 and supplementary note ) . these predictive power levels were obtained using the full history models , that is , the instantaneous spiking rate of each neuron conditioned on both intrinsic and ensemble histories . we further compared the predictive power of each of these two components separately ( fig . the predictive power obtained exclusively from the ensemble histories was larger than the predictive power from intrinsic histories . in addition , it is possible that some of the intrinsic history effects , beyond refractory and recovery period effects , could also reflect network dynamics . thus , these results indicate that the collective ensemble activity , rather than refractory and recovery dynamics , was the main source of the observed predictive power of spiking histories . the predictive power of these m1 spiking history models in monkeys mla and mcl was also higher than the predictive power of relevant kinematic covariates , such as hand position and velocity ( supplementary fig . 3 and supplementary note ) . we used pathlet models21 to estimate spiking probabilities conditioned on kinematics ; that is , the instantaneous spiking rate was expressed as function of hand position and velocity trajectory in the time interval ( 200 ms , 300 ms ) with respect to spiking time . furthermore , the predictive power values of spiking history and pathlet models were not strictly additive ( supplementary fig . this result indicates that there was some level of redundancy between the information conveyed by these two models about single - neuron spiking : some of the spiking activity that was predicted by the spiking history models could also be predicted by the examined kinematics . beyond examining the predictive power of ensemble spiking histories , we were also interested in the predictive power of instantaneous collective states in the recorded ensemble , that is , simultaneous states at either 1- or 10-ms temporal resolution . instantaneously correlated states could result from common inputs and/or from synchronization patterns arising from the neuronal network s own dynamics . strong instantaneous collective states7 could still be consistent with the weak pair - wise instantaneous interactions that we observed in our datasets ( fig . we estimated the joint probability distribution of these instantaneous collective states by fitting maximum entropy22,23 point process models constrained by empirical mean spiking rates and zero time - lag pair - wise correlations . that is , these maximum entropy distribution models were consistent with the observed mean spiking rates and zero time - lag pair - wise correlations , but made no additional assumptions ( see online methods ) . this pair - wise maximum entropy probability model , also known as the ising model in statistical mechanics24,25 , has been shown to capture most of the instantaneous interdependency structure in ex vivo / in vitro neuronal ensemble preparations7,8,10 . besides satisfying the pair - wise correlation structure in the examined datasets , the estimated models also accounted well for the distribution of multiple - neuron spike coincidences , as seen in the comparison of the empirical distribution and the distribution generated from the model via gibbs sampling26 ( fig . , the predictive power of this pair - wise correlation structure was considerably poorer than the predictive power of intrinsic and ensemble spiking histories , as captured by the models ( fig . - wise maximum entropy models corresponded to 0.10  0.17 , while the spiking history models yielded values of 0.50  0.40 ( mean  s.d . ) . the predictive power of spiking history models were still higher when we considered instantaneous states at a coarser 10-ms temporal resolution ( fig . 4 ) , which resulted in predictive power values of 0.26  0.27 ( mean  s.d . ) . we performed the same predictive power analysis for areas pmv and 5d and compared predictions based on ensemble histories recorded from the same or a different cortical area in the connected pairs m1pmv and m15d ( fig . predictive power , based on both local and inter - areal spiking histories , corresponding to target neurons in pmv ( mco ; fig . 5c ) had medians of 0.27 and 0.32 and ranges of 0.110.86 and 0.130.62 , respectively . predictive power corresponding to target neurons in m1 , based on both local and inter - areal spiking histories , had medians of 0.50 and 0.40 and ranges of 0.150.93 and 0.110.83 for monkeys mco and mab , respectively ( fig . the inter - areal predictive power could also be substantial in both directions , m1pmv and m15d ( fig . medians were 0.15 , 0.13 , 0.13 and 0.09 and ranges were 0.050.73 , 0.060.75 , 0.060.37 and 0.050.37 for m1pmv , pmvm1 , m15d and 5dm1 , respectively . furthermore , the predictive power from one area to target neurons in the other , especially from m1pmv , could be larger than the predictive power of local ensembles ( fig . we speculated that the higher predictive power of m1pmv in some of these cases could simply reflect the larger number of recorded neurons in m1 , rather than the relative strength of feedforward and feedback connections between these two areas . a supplementary analysis using ensembles of equal size suggested that this was the case ( supplementary fig . 5 and supplementary note ) . this analysis also indicated that intra - areal predictive power was slightly higher than inter - areal predictive power when ensembles of equal size were considered . as shown before for subjects mla and mcl , the predictive power of spiking history models for areas m1 , 5d and pmv also tended to be higher than the predictive power of pathlet models ( supplementary fig . the medians and ranges of predictive power values for pathlet models corresponded to 0.25 and 0.21 , and 0.010.77 and 0.010.75 for mco and mab , respectively . the predictive power of spiking history models varied broadly across the 1,187 neurons in the 12 datasets . we examined whether this variation could be easily explained by simple features of the predicted spiking activity . the predictive power of history models did not appear to depend , at least in a simple manner , on the mean spiking rate or on the level of irregularity of the predicted spiking activity ( fig . assessments of spike train irregularity were based on the coefficient of variation3 of the interspike time intervals . we also assessed how much information was involved in the prediction of single neuron spiking . in principle , the same level of prediction accuracy could be achieved for two different target neurons while involving different amounts of information . we examined this possibility by asking how much information about a neuron spiking or not at any given time was gained by modeling the intrinsic and ensemble history effects compared with knowing only the mean spiking rate of the process . this information gain can be estimated as the normalized difference between the log - likelihood under the estimated history model and the log - likelihood under a homogenous poisson process with given mean spiking rate ( see online methods ) . our analysis confirmed that similar predictive power values ( in the prediction of different target neurons ) could involve a wide range of information rates ( fig . on the basis of the estimated conditional spiking probabilities , we used a standard tool , the roc curve analysis19 , to assess the predictive power of spiking history models ( see online methods ) . intra - areal ensemble history  based predictions in m1 resulted in high true - positive rates while maintaining low numbers of false positives . for example , it was possible to correctly predict 80% of spikes in neuron 34a ( participant hs3 ) with a false - positive rate of less than 5% ( fig . similarly , an 80% true - positive rate was obtained for neuron 16a ( monkey mla ) with a less than 10% false - positive rate ( fig . this predictability was substantial despite the weak pair - wise correlations that we observed among all neuronal pairs in the recorded datasets for hs3 and mla ( fig . 2c , we obtained a more comprehensive assessment of predictive performance by computing the area under the roc curve . the area under the curve ( auc ) is a global summary statistic ; that is , it depends on both the true- and false - positive rates and on all of the possible thresholds on the spiking probability . the auc gives the probability that , when two samples are randomly drawn from the data ( one containing a spike , the other not ) , the conditional intensity model will assign a higher probability ( i.e. a higher instantaneous spiking rate ) to the sample with a spike20 . it therefore provides an assessment of the discriminatory power for predictive variables under a given model . it approaches 0.5 for a chance level predictor , that is , a predictor having false- and true - positive rates along the diagonal , and it equals 1 for a perfect predictor . for example , the auc for neuron 34a was 0.95 ( fig . we further corrected the auc by subtracting chance - level predictions estimated from the actual data and scaled it to obtain a quantity that ranged from 0 ( no predictive power ) to 1 ( perfect prediction ) . we refer to this corrected and normalized auc measure as the predictive power ( see online methods ) . more than 50% of examined neurons in monkey mla showed a predictive power higher than 0.5 ( median = 0.55 , range = 0.10.89 ) ( fig . 3 ) . prediction in subjects mcl and hs3 showed similar performances , with medians of 0.45 and 0.49 and ranges of 0.180.91 and 0.230.97 , respectively . predictive power was lower for participant hs1 ( median = 0.30 , range = 0140.69 ) . we conjectured that this was probably a result of the smaller number of neurons ( n = 21 , n = 22 ) that were recorded and used for prediction in this participant . to test this conjecture , we computed the predictive power on the basis of several smaller ensemble subsets ( n = 22 neurons ) randomly chosen out of the 110 neurons from hs3 s dataset ( session 2 ) . data from hs3 provided a good reference , as neurons were recorded from m1 under the same task condition as for hs1 . the distribution of predictive power values , based on this collection of smaller ensemble subsets , was similar to the one obtained for hs1 , supporting our hypothesis ( supplementary fig . 2 and supplementary note ) . these predictive power levels were obtained using the full history models , that is , the instantaneous spiking rate of each neuron conditioned on both intrinsic and ensemble histories . we further compared the predictive power of each of these two components separately ( fig . the predictive power obtained exclusively from the ensemble histories was larger than the predictive power from intrinsic histories . in addition , it is possible that some of the intrinsic history effects , beyond refractory and recovery period effects , could also reflect network dynamics . thus , these results indicate that the collective ensemble activity , rather than refractory and recovery dynamics , was the main source of the observed predictive power of spiking histories . the predictive power of these m1 spiking history models in monkeys mla and mcl was also higher than the predictive power of relevant kinematic covariates , such as hand position and velocity ( supplementary fig . 3 and supplementary note ) . we used pathlet models21 to estimate spiking probabilities conditioned on kinematics ; that is , the instantaneous spiking rate was expressed as function of hand position and velocity trajectory in the time interval ( 200 ms , 300 ms ) with respect to spiking time . furthermore , the predictive power values of spiking history and pathlet models were not strictly additive ( supplementary fig . this result indicates that there was some level of redundancy between the information conveyed by these two models about single - neuron spiking : some of the spiking activity that was predicted by the spiking history models could also be predicted by the examined kinematics . beyond examining the predictive power of ensemble spiking histories , we were also interested in the predictive power of instantaneous collective states in the recorded ensemble , that is , simultaneous states at either 1- or 10-ms temporal resolution . instantaneously correlated states could result from common inputs and/or from synchronization patterns arising from the neuronal network s own dynamics . strong instantaneous collective states7 could still be consistent with the weak pair - wise instantaneous interactions that we observed in our datasets ( fig . we estimated the joint probability distribution of these instantaneous collective states by fitting maximum entropy22,23 point process models constrained by empirical mean spiking rates and zero time - lag pair - wise correlations . that is , these maximum entropy distribution models were consistent with the observed mean spiking rates and zero time - lag pair - wise correlations , but made no additional assumptions ( see online methods ) . this pair - wise maximum entropy probability model , also known as the ising model in statistical mechanics24,25 , has been shown to capture most of the instantaneous interdependency structure in ex vivo / in vitro neuronal ensemble preparations7,8,10 . besides satisfying the pair - wise correlation structure in the examined datasets , the estimated models also accounted well for the distribution of multiple - neuron spike coincidences , as seen in the comparison of the empirical distribution and the distribution generated from the model via gibbs sampling26 ( fig . , the predictive power of this pair - wise correlation structure was considerably poorer than the predictive power of intrinsic and ensemble spiking histories , as captured by the models ( fig . 4 ) . across the four subjects , the predictive power values of pair - wise maximum entropy models corresponded to 0.10  0.17 , while the spiking history models yielded values of 0.50  0.40 ( mean  s.d . ) . the predictive power of spiking history models were still higher when we considered instantaneous states at a coarser 10-ms temporal resolution ( fig . 4 ) , which resulted in predictive power values of 0.26  0.27 ( mean  s.d . ) . we performed the same predictive power analysis for areas pmv and 5d and compared predictions based on ensemble histories recorded from the same or a different cortical area in the connected pairs m1pmv and m15d ( fig . predictive power , based on both local and inter - areal spiking histories , corresponding to target neurons in pmv ( mco ; fig . 5c ) had medians of 0.27 and 0.32 and ranges of 0.110.86 and 0.130.62 , respectively . predictive power corresponding to target neurons in m1 , based on both local and inter - areal spiking histories , had medians of 0.50 and 0.40 and ranges of 0.150.93 and 0.110.83 for monkeys mco and mab , respectively ( fig . the inter - areal predictive power could also be substantial in both directions , m1pmv and m15d ( fig . medians were 0.15 , 0.13 , 0.13 and 0.09 and ranges were 0.050.73 , 0.060.75 , 0.060.37 and 0.050.37 for m1pmv , pmvm1 , m15d and 5dm1 , respectively . furthermore , the predictive power from one area to target neurons in the other , especially from m1pmv , could be larger than the predictive power of local ensembles ( fig . we speculated that the higher predictive power of m1pmv in some of these cases could simply reflect the larger number of recorded neurons in m1 , rather than the relative strength of feedforward and feedback connections between these two areas . a supplementary analysis using ensembles of equal size suggested that this was the case ( supplementary fig . 5 and supplementary note ) . this analysis also indicated that intra - areal predictive power was slightly higher than inter - areal predictive power when ensembles of equal size were considered . as shown before for subjects mla and mcl , the predictive power of spiking history models for areas m1 , 5d and pmv also tended to be higher than the predictive power of pathlet models ( supplementary fig . the medians and ranges of predictive power values for pathlet models corresponded to 0.25 and 0.21 , and 0.010.77 and 0.010.75 for mco and mab , respectively . the predictive power of spiking history models varied broadly across the 1,187 neurons in the 12 datasets . we examined whether this variation could be easily explained by simple features of the predicted spiking activity . the predictive power of history models did not appear to depend , at least in a simple manner , on the mean spiking rate or on the level of irregularity of the predicted spiking activity ( fig . assessments of spike train irregularity were based on the coefficient of variation3 of the interspike time intervals . we also assessed how much information was involved in the prediction of single neuron spiking . in principle , the same level of prediction accuracy could be achieved for two different target neurons while involving different amounts of information . we examined this possibility by asking how much information about a neuron spiking or not at any given time was gained by modeling the intrinsic and ensemble history effects compared with knowing only the mean spiking rate of the process . this information gain can be estimated as the normalized difference between the log - likelihood under the estimated history model and the log - likelihood under a homogenous poisson process with given mean spiking rate ( see online methods ) . our analysis confirmed that similar predictive power values ( in the prediction of different target neurons ) could involve a wide range of information rates ( fig . our findings reveal that single - neuron spiking in the cortex of humans and monkeys performing sensorimotor tasks can be substantially predicted by the spiking history of small , randomly sampled neuronal ensembles . the fact that this predictability was based on small neuronal ensembles , randomly sampled out of millions of neurons in the cortex , suggests that strongly coordinated activity underlies the generation of single - neuron spikes . cortical neurons are embedded in large , sparsely connected recurrent networks in which the high number of synaptic inputs and high - conductance states typically induce weak coupling between randomly sampled neuron pairs . not surprisingly , revealing and understanding these large scale and dynamic interactions has been challenging . on the basis of only the weak pair - wise correlations observed amongst cortical neurons in our datasets , one would have underestimated the strength of the statistical interdependencies induced by the collective dynamics . furthermore , we believe our estimates of predictive power for these small neuronal ensembles should be taken as lower - bound values . for example , even though we were careful to include only trials or time segments in which the hand was moving or when the human participants were controlling a computer cursor , it is possible that the network s functional connectivity was nonstationary within and across trials . also , we avoided more complex , and thus potentially more predictive , point process history models for computational tractability . the fact that spiking was better predicted by the ensemble s history than by the ensemble s simultaneous collective state , estimated via pair - wise maximum entropy models , emphasizes the temporal dynamics leading to spiking . this finding , however , should not be taken as a limitation of pair - wise maximum entropy models . it is possible that multiple time - lag pair - wise correlation maximum entropy models11,27 might capture most of the history effects detected in our data and therefore provide a simpler , minimal model . our goal here was not , however , to provide such a minimal model for the temporal and instantaneous collective dynamics , but to determine the existence and strength of such dynamics in the cortex of humans and behaving monkeys . this is , nevertheless , an important issue that should be addressed . in previous studies , we and others have shown that the spiking activity of small , simultaneously and randomly sampled neuronal ensembles in motor cortex can be used to predict ( decode ) subsequent complex behavioral variables such as arm kinematics2832 . here , we found that the activity of these same neuronal ensembles can also be used to predict subsequent single - neuron spiking with substantial accuracy , implying the presence of strong collective dynamics in sensorimotor cortex . although a comprehensive analysis of this problem is beyond the scope of this study , our results indicate that these collective dynamics do not simply reflect background coherent states that are completely unrelated to behavior and , conversely , do not simply reflect  trivial common inputs , such as those usually considered in studies involving stimulus - driven activities of early sensory neurons . regarding the background coherent states , our data indicate that the predictive power of models based only on ensemble history and of models based only on kinematics are not strictly additive ( supplementary fig . , there is some level of redundancy between the information ( about single - neuron spiking ) conveyed by ensemble spiking histories and information conveyed by kinematics variables . therefore , the detected collective dynamics can not simply reflect background coherent states that are entirely unrelated to behavior . nonetheless , it is still possible , if not likely , that part of the detected collective dynamics may reflect ongoing internal processes that are not related to behavioral variables . regarding common inputs , we note that kinematics and other features of voluntary movements are controlled in large part ( either directly , or indirectly via spinal and muscle activations ) by the coordinated activity of sensorimotor cortical neurons . it is , therefore , not surprising that these behavioral variables can also achieve substantial explanatory power for spiking activity . even so , the fact that kinematics ( pathlet ) models were less predictive than ensemble history models suggests that the latter carried extra information about single - neuron spiking ( supplementary fig . the detected collective dynamics are unlikely to be explained as the simple reflection of trivial common inputs . in sum , we believe that these strong collective dynamics reflect the intra and inter - areal coordinated activity of neuronal ensembles distributed in the many different cortical and subcortical areas that participate in voluntary control of movement . we hypothesize that the detected collective dynamics and ensemble influences on spiking activity reflect information transfer and computation in cortical networks . collective dynamics and functional connectivity , as captured by connectivity matrices derived from ensemble history models , as well as predictive power levels , should vary as information transfer and computation change during behaviors that engage cortical areas . on the basis of current computational theories of motor control14,15 , one could predict , for example , that m15d spiking predictability will manifest primarily during the initial phase of reaching movements , whereas m1pmv spiking predictability will peak during the hand - shaping phase of object grasping . a related and more general inquiry will be to examine the relationship of collective dynamics at this microscopic spatial scale to neural activities reflected in meso- and macroscopic scale signals , such as local field potentials and electrocorticograms . our results also have implications for neural decoding theory and intracortical neural interfaces for motor prostheses . collective dynamics add redundancy and , therefore , error - correcting properties to neural codes7 . in addition , these dynamics might also account for variability of neural responses5 , which is otherwise usually attributed to noise . therefore , it seems that ensemble history effects should be taken into consideration when decoding kinematics ( or other variables ) from the spiking activity of neuronal ensembles . one would predict that decoding algorithms that take into account ensemble spiking histories will outperform algorithms that treat spiking activity of different neurons as , conditioned on decoded variables , independent processes . our findings suggest the presence of strong collective dynamics that are fundamental to the experimental endeavor of determining coordinated spiking in cortical networks . networks responsible for specialized cortical function are likely to be contained in the spiking patterns of millions of neurons distributed across multiple cortical areas . current and developing recording technologies measure the spiking activity of hundreds to thousands of neurons , a very small fraction of these networks . without strong collective dynamics ( that is , if neurons in small randomly sampled ensembles behaved seemingly independently ) , there would be little hope of determining how the coordinated propagation of action potentials in large - scale recurrent networks leads to computation and information processing . we believe , therefore , that the existence of these collective dynamics offers a basis for understanding cognition and adaptive behavior at the level of coordinated spiking in cortical networks . methods and any associated references are available in the online version of the paper at http://www.nature.com / natureneuroscience/.