A Granger Causality Measure for Point Process Models of Ensemble Neural Spiking Activity

The ability to identify directional interactions that occur among multiple neurons in the brain is crucial to an understanding of how groups of neurons cooperate in order to generate specific brain functions. However, an optimal method of assessing these interactions has not been established. Granger causality has proven to be an effective method for the analysis of the directional interactions between multiple sets of continuous-valued data, but cannot be applied to neural spike train recordings due to their discrete nature. This paper proposes a point process framework that enables Granger causality to be applied to point process data such as neural spike trains. The proposed framework uses the point process likelihood function to relate a neuron's spiking probability to possible covariates, such as its own spiking history and the concurrent activity of simultaneously recorded neurons. Granger causality is assessed based on the relative reduction of the point process likelihood of one neuron obtained excluding one of its covariates compared to the likelihood obtained using all of its covariates. The method was tested on simulated data, and then applied to neural activity recorded from the primary motor cortex (MI) of a Felis catus subject. The interactions present in the simulated data were predicted with a high degree of accuracy, and when applied to the real neural data, the proposed method identified causal relationships between many of the recorded neurons. This paper proposes a novel method that successfully applies Granger causality to point process data, and has the potential to provide unique physiological insights when applied to neural spike trains.


Introduction
Neurons in the brain are known to exert measurable, directional influences on the firing activities of surrounding neurons, and a detailed analysis of these interactions improves our understanding of how the brain performs specific functions [1]. Attempts to identify associations between neurons, such as the cross-correlogram [2], joint peri-stimulus time histogram [3], smoothed ratio of spiking activity [4], and gravitational clustering [5], have been useful in the past. However, these methods provide little insight into the directional nature of the interactions that they detect, are less reliable in their detection of inhibitory interactions, and usually do not consider the point process nature of neural spike train data. Occasionally they may also give a misleading picture of the relationships between neurons if the detected associations are caused by common inputs or mediated by other neurons [6].
Granger causality has proven to be an effective method for the investigation of directional relationships between continuousvalued signals in many applications [7][8][9][10][11]. The basic idea of causality between signals was introduced by Wiener [12] but was too general to be implemented. Granger formalized this idea in order to enable practical implementation based on the multivariate autoregressive (MVAR) models [7]: if past values of y contain information that helps predict x above and beyond the information contained in past values of x alone, then y is said to Granger-cause (or g-cause) x. Its mathematical formulation is based on the MVAR modeling of processes. However, it is difficult to apply this method directly to spike train data, since they can not be described by the MVAR model, and standard distance measures such as the mean squared error (MSE) are not designed for spike train data. Recently, several methods have been developed to apply Granger causality to spike train data [13][14][15][16][17][18][19]. Attempts at transforming neural spike trains into continuousvalued data by convolving spike trains with either a smooth kernel [13] or a lowpass filter [14,15] have been proposed, but they introduced unwanted distortion of the point process characteristics of spike train data. Granger causality analysis based on an MVARnonlinear-Poisson model has been proposed [16]; however, this approach lacks an explanation of the physical meaning of the model that is being applied. A method called transfer entropy using mutual information has also been proposed [17,18], and it is sensitive to nonlinear signal properties, but unfortunately its application is restricted to bivariate cases. A nonparametric method based on spectral matrix factorization has been proposed [19]; however, it required the second-order stationarity of spike train data.
To address these issues, this paper proposes a point process framework for assessing Granger causality between multiple neurons. The spiking activity of each neuron is simultaneously affected by multiple covariates such as its own spiking history and the concurrent ensemble activity of other neurons. The effect of these factors on a neuron's spiking activity is characterized by a statistical framework based on the point process likelihood function, which relates the neuron's spiking probability to the covariates [20,21]. Using the point process likelihood function, Granger causality between neurons is assessed based on the likelihood ratio statistic. That is, Granger causality from neuron j to neuron i is measured based on the relative reduction of the point process likelihood of neuron i obtained by excluding the covariates corresponding to the effect of neuron j compared to the likelihood obtained using all the covariates. If the likelihood ratio is less than one, we say that there is a causal influence from neuron j to i, and if the ratio is one, we say that there is no causal influence. In continuous-valued cases, the Granger causality measure based on the MVAR prediction error was shown to be the likelihood ratio test statistic if the prediction error is assumed to be Gaussian [22]. In addition, the point process likelihood ratio statistic enables us to perform statistical hypothesis testing to investigate the significant causal interactions between neurons, since it asymptotically follows a chi-squared distribution when the conditional intensity function (CIF) of the point process is modeled by the generalized linear model (GLM) [23]. When performing a set of statistical inferences simultaneously to detect statistically significant causal interactions among all possible interactions, multiple hypothesis testing problems where the null hypothesis is more likely to be incorrectly rejected should be considered. The present study uses the false discovery rate (FDR) correction to control the expected proportion of incorrectly rejected null hypotheses [24].
The proposed framework was used in an attempt to identify the causal relationships between simulated spike train data, and accurately estimated the underlying causal networks presented in the simulations. It was also applied to real neural data recorded from the cat primary motor cortex (MI) in order to assess the causal relationships that occur between multiple simultaneously recorded neurons during performance of a movement task.

Ethics Statement
The experiments that were performed for the collection of real neural spiking data were approved by the Animal Ethics Committee of the University of Western Australia, and the National Health and Medical Research Council of Australia (NH&MRC) guidelines for the use of animals in experiments were followed throughout.

Summary of the Proposed Method
Statistical analysis of the potential causal relationships between neurons was performed based on a point process likelihood function. The likelihood function related a neuron's spiking probability to possible covariates, such as its own spiking history and the concurrent activity of all simultaneously recorded neurons. The causal relationships between associated neurons were assessed based on the point process likelihood ratio, which represents the extent to which the likelihood of one neuron is reduced by the exclusion of one of its covariates, compared with the likelihood if all of the available covariates are used. The Granger causality measure based on the point process likelihood ratio also enabled us to detect significant causal relationship through a hypothesis testing based on the likelihood ratio statistic.

A Granger Causality Measure for Point Process Models
A point process is a time series of discrete events that occur in continuous time [25]. The discrete, all-or-nothing nature of a sequence of action potentials together with their stochastic structure suggests that neural spike trains may be regarded as point processes [20,[26][27][28]. Given an observation interval (0,T, ƒT be a set of J i spike times point process observations for i~1,:::,Q recorded neurons. Let N i (t) denote the sample path that counts the number of spikes of neuron i in the time interval (0,t for t [ (0,T. A point process model of a spike train for neuron i can be completely characterized by its CIF, l i (tjH i (t)), defined as where H i (t) denotes the spiking history of all the neurons in the ensemble up to time t for neuron i [25]. In this work, H i (t) is defined in the interval ½t{M i W ,t), which is divided into M i nonoverlapping rectangular windows of duration W ; We denote the spike count of neuron q in a time window of length W covering the time interval ½t{mW ,t{(m{1)W ) as R q,m (t) for q~1,:::,Q and m~1,:::,M i . The CIF, l i (tjH i (t)), of (1) represents the firing rate of neuron i at time t, so it quantifies the probability that neuron i fires a spike at time t given its covariates H i (t). Each neuron has a different H i (t), since each has a history dependency of different length, M i W . The probability that neuron i fires a single spike in a small interval ½t,tzD) can be approximated as l i (tjH i (t))D.
To model the effect of its own and ensemble's spiking histories on the current spiking activity of a neuron, a GLM framework is often used to model the CIF. In the GLM framework, the logarithm of the CIF is modeled as a linear combination of the functions of the covariates that describe the neural activity dependencies [20,21]. Thus, the logarithm of the CIF is expressed as where c i,0 relates to a background level of activity of neuron i, and

Author Summary
Recent advances in multiple-electrode recording have made it possible to record the activities of multiple neurons simultaneously. This provides an opportunity to study how groups of neurons form functional ensembles as different brain areas perform their various functions. However, most of the methods that attempt to identify associations between neurons provide little insight into the directional nature of the interactions that they detect. Recently, Granger causality has proven to be an efficient method to infer causal relationships between sets of continuous-valued data, but cannot be directly applied to point process data such as neural spike trains. Here, we propose a novel and successful attempt to expand the application of Granger causality to point process data. The proposed method performed well with simulated data, and was then applied to real experimental data recorded from sets of simultaneously recorded neurons from the primary motor cortex. The results of the real data analysis suggest that the proposed method has the potential to provide unique neurophysiological insights about network properties in the cortex that have not been possible with other contemporary methods of functional interaction detection.
c i,q,m represents the effect of ensemble spiking history R q,m (t) on the firing probability of neuron i. The parameter vector c i is given as c i~f c i,0 ,c i,1,1 ,:::,c i,q,m ,:::,c i,Q,M i g, ð3Þ which represents the dependency of neuron i on the spiking history of all neurons in the ensemble. Especially, the parameters fc i,q,m g M i m~1 represent the dependency of neuron i on the spiking history of neuron q for i,q~1,:::,Q. The model for the CIF of (2) is not a fixed form, but can change depending on its covariates and its relationship to them.
A point process likelihood function was used to fit the parametric CIF and analyze Granger causality between neurons since it is a primary tool used in constructing statistical models and has several optimality properties [29]. Here, we used a discrete time representation of the point process likelihood function in order to simplify ensuing calculations. To obtain this representation, we partitioned the observation interval (0,T into K subintervals (t k{1 ,t k K k~1 each of length D~TK {1 where K is a large integer. Usually, K is chosen to make D as 1 ms. We denote the continuous time variables defined above as the discrete time versions such as N i ½k for N i (t k ), H i ½k for H i (t k ), R q,m ½k for R q,m (t k ) and so forth. Since we chose a large value for K, there is at most one spike per subinterval, that is, DN i ½kÑ i ½k{N i ½k{1 takes on the value 0 if there is no spike in (t k{1 ,t k or 1 if there is a spike. The parametric form of the CIF of (2) for neuron i is represented as l i (t k jc i ,H i ½k).
Given the ensemble spiking activity in (0,T, the likelihood function of the spike train of neuron i is given as in [20] using its CIF by where the term o(D J i ) relates the probability that neuron i includes two or more spikes in any subinterval (t k{1 ,t k . Based on the likelihood function of (4), a point process framework for assessing the causal relationships between neurons is proposed. A potential causal relationship from neuron j to neuron i is assessed by calculating the relative reduction in the likelihood of producing a particular set of spike trains of neuron i if the spiking history of neuron j is excluded, compared with the likelihood if all of the available covariates are used. The log-likelihood ratio, C ij , is given by where the likelihood L i (c j i ) is obtained using a new CIF, l j i , which excludes the effect of neuron j from H i ½k, given as The parameter vector c j i is obtained by re-optimizing the parametric likelihood model after excluding fc i,j,m g M i m~1 from c i in order to remove the effect of neuron j on neuron i, and H j i ½k is obtained by leaving out fR j,m ½kg M i m~1 from H i ½k. Since the likelihood L i (c i ) is always greater than or equal to the likelihood L i (c j i ), the log-likelihood ratio C ij is always less than or equal to 0. If the spiking activity of neuron j has a causal influence on that of neuron i in the Granger sense, the likelihood L i (c i ) that is calculated using all the covariates of neuron i is greater than the likelihood L i (c j i ) that is calculated using the same covariates, save for the history of neuron j, which is excluded. Excitatory and inhibitory influences of neuron j on neuron i can be distinguished by the sign of P M i m~1 c i,j,m that represents an averaged influence of the spiking history of neuron j on neuron i. The equality holds when neuron j has no influence on neuron i. Thus, the Granger causality measure from neuron j to neuron i is proposed as which provides an indication of the extent to which the spiking history of neuron j affects the spike train data of neuron i. A positive result is indicative of neuron j having an excitatory effect upon neuron i, a negative result indicates an inhibitory effect, and zero indicates that no interactions are detected. Finally, a Q|Q Granger causality matrix can be produced, W, whose (i,j)th element is w ij , and represents the extent to which neuron j has either an excitatory or inhibitory influence on neuron i for i,j~1,:::,Q.

Significance Test
The Granger causality matrix W represents the relative strength of estimated causal interactions between neurons, but does not provide any insight into which of these interactions are statistically significant. To address this issue, a hypothesis testing based on the likelihood ratio test statistic is performed to evaluate the statistical significance of the estimated causal interactions of W. For this, the goodness-of-fit (GOF) statistics are applied as follows [23,29]. Let us denote the deviance obtained using the model parameter c j i as D 0 and the deviance obtained using the model parameter c i as D 1 ; The deviance is obtained by comparing the estimated model with a more general model that has a parameter for every observation so that the data fits exactly, which is called a full model [25,30]. Its expression is 22 times the log-likelihood ratio of the estimated model to the full model, which is mathematically expressed by where c and c max are the parameters for the estimated and the full models, respectively. In the GLM framework the deviance is used to compare two models, which are nested like the above case, since a model of c j i is a special case of the more general model of c i . Consider the null hypothesis which corresponds to the model of (6). An alternative hypothesis is which corresponds to the model of (2). We can test H 0 against H 1 using the difference of the deviance statistic as the test statistic, which is given by Thus, the deviance difference between two models is equivalent to 22 times log-likelihood ratio given by (5). If both models describe the data well, then the deviance difference may be asymptotically described as DD*x 2 M i where M i is equal to the difference in dimensionality of the two models [23,29]. If the value of DD is consistent with the x 2 M i distribution, the hypothesis H 0 is accepted since it is simpler. This result indicates that the past values of neuron j contain no significant information that would assist in predicting the activity of neuron i. Thus, neuron j has no causal influence on neuron i. On the contrary, if the value of DD is in the critical region, i.e., greater than the upper tail 100|(1{a)% of the x 2 M i distribution where a determines false positive rates, then H 0 may be rejected in favor of H 1 since the model of (2) describes the data with significantly more accuracy. This indicates that past spike times of neuron j contain information that improves the ability to predict the activity of neuron i. Thus the activity of neuron j g-causes the activity of neuron i.
In any attempt to identify the causal relationships between multiple neurons simultaneously, the total number of the possible causal interactions to be investigated is usually large. Thus, the use of common statistical thresholds cited above to assess the causal interactions would lead to an unacceptably large number of false causal interactions (false positives) where the null hypothesis is incorrectly rejected [31]. The multiple comparison problem could potentially be addressed by the use of stricter statistical thresholds, which would result in a reduction in the proportion of the falsely rejected null hypotheses. However, stricter thresholds would also reduce the probability that true causal interactions between neurons were identified. The present study uses a multiplehypothesis testing error measure called the FDR to address the multiple comparisons problem. The FDR controls the expected proportion of false positive findings among all the rejected null hypotheses [24]. In situations where the number of hypothesis tests is large, other approaches that attempt to control the familywise error rate (FWER), which is the probability of making one or more false discoveries among all the hypotheses, can be too strict and decrease the power. Thus, the FDR is a less conservative, but more powerful, quantity to control for multiple comparisons than the FWER at a cost of increasing the likelihood of obtaining false positive findings [32].
Combining the multiple hypothesis testing results with the sign of P Mi m~1 c i,j,m , we detect the inhibitory, excitatory, and non-causal interactions, which are denoted as the blue, red, and green colors, respectively. Thus, a Q|Q causal connectivity matrix Y whose (i,j)th element corresponds to one of three interactions is constructed. In this paper, the connectivity matrix Y was obtained by controlling the FDR as 0.05.

Simulation
In order to evaluate the proposed framework's ability to identify Granger causality for ensemble spiking activity, we analyzed synthetically generated spike train data. Simulated spike train data were synthetically generated based on the nine-neuron network of Figure 1. The firing probability of each neuron was dependent on its own spiking history and the concurrent ensemble activity through the inhibitory and the excitatory interactions of Figure 1. The inhibitory and the excitatory interactions were represented as black and white circles, respectively. Each neuron was influenced by other neurons through two inhibitory interactions including its own self-inhibition and through one or two excitatory interactions. The overall network of Figure 1A consisted of three sub-networks All neurons had the same spontaneous firing rate (18 Hz). Spike trains for neuron i were generated using a commonly used procedure as follows [33]: A random number r, uniformly distributed between 0 and 1, is generated at every interval; if rƒl i (t k jc i ,H i ½k)D, a spike is presumed to have occurred in (t k ,t k zD; otherwise, no spike is generated. The time resolution D was set to 1 ms. An absolute refractory period of 1 ms was enforced to prevent neurons from firing a spike in adjacent time steps. Based on the experimental settings cited above, we generated 100,000 samples for each neuron, and the total number of spikes for each neuron ranged from 2176 through 2911. Examples of generated neural spike trains during the first 5 sec (5,000 samples) are illustrated in Figure 2. It can be seen that neurons generally fire less (or more) spikes after other neurons with inhibitory (or excitatory) influence on them fire spikes. However, it is hard to estimate the underlying causal network between neurons from this plot.
In order to select a model for each neuron we fit several models with different history durations M i W to each spike train data and then identified the best approximating model from among a set of candidates using Akaike's information criterion (AIC) [34,35]. Using this criterion, an optimum model order for each neuron was selected. The spike counting window length W was set to 2 ms.  For neurons 1, 3, 4, 5, 8, and 9, which were influenced by other neurons through relatively long interactions, the selected GLM spike order M i was 3, which indicates a 6 ms history duration, and for neurons 2, 6, and 7 influenced by other neurons within same sub-network only through short interactions, the selected GLM order M i was 2, which corresponds to a 4 ms history duration.
Based on the estimated model, two kinds of causality maps were obtained using the proposed method. Firstly, the Granger causality map W, which is illustrated in Figure 3A, represents the relative strength of the causal interaction between neurons. It represents the extent to which a trigger neuron has a causal impact on a target compared to other interconnections, but provides little insight into which causal impact is statistically significant. In order to make up for W, the causal connectivity map Y was obtained through the hypothesis testing when we controlled the FDR as 0.05. This is shown in Figure 3B. The red, blue, and green colors denote the presence of excitatory, inhibitory, or no interactions from trigger neuron to target, respectively. The estimated pattern of Y matches the actual network of Figure 1 exactly. This causality map does not show a connection between neurons that do not have direct interactions, even though they have indirect interactions.
The FDR procedure was used as a solution for the multiple comparisons problem when considering a set of statistical inferences simultaneously. When controlling the FDR at a specific significance level k, we expect that on average there will be kR false positives amongst R detected significant interactions. In order to verify that the FDR is actually being controlled at the significance levels that we are claiming in the present study, the Monte-Carlo (MC) simulations were conducted by varying both the number of causal interactions between the virtual neurons, and the signal-to-noise ratio (SNR) of the simulated spikes. These MC simulations show how effectively the FDR is being controlled under different experimental conditions. Firstly, we conducted a series of the MC simulations by changing the number of causal interactions from 8 to 64. Data was synthetically generated to resemble four different kinds of networks (seen in Figure 4), each having a different incidence of interaction between neurons. Fifty data sets were generated for each network condition, while all other experimental parameters remained the same. The dashed circle in Figure 4A and B represents a neuron whose firing activity does not depend on the spiking history, and thus follows a homogeneous Poisson process, i.e., c (0) i,i, : = 0. Networks of Figure 4A to D consist of 8, 16, 32, and 64 interactions (including selfinteractions), respectively. The observed FDR is calculated by averaging the ratio of the number of false positives to the number of detected significant interactions over 50 simulations, and it is illustrated in Figure 5A for significance levels of 0.01, 0.05 and 0.1. The FDR was generally controlled at the significance level that we were attempting to control except for the 8-interaction case with less false positives than the number expected at that significance level. We then performed another MC simulation by changing the SNR. Noisy neural spike trains were generated using the CIF of (2) in the following: We added a Gaussian noise to the logarithmic CIF, i.e., the right-hand side of (2), and then generated spike trains using the perturbed CIF. The noise changed the background level of firing rate over time. The SNR is defined as the ratio between the unperturbed logarithmic CIF and the perturbation itself. Fifty data sets of noisy spike trains were synthetically generated based on the nine-neuron network of Figure 1 with different levels of noise, which led to about 0, 10, 20, 30, and 40 dB SNRs, respectively. All other experimental conditions are same to the previous case. Figure 5B illustrates the simulation results obtained for significance levels of 0.01, 0.05, and 0.1. When the SNR is approximately 0 dB, more false positive events were detected than what was expected at the specified significance level, but in most cases the observed FDR was no different from the theoretical FDR. In summary, unless the perturbation level is similar to or higher than the level of the logarithmic CIF of (2) that is modulated by the intrinsic dynamics of the neurons, the FDR is effectively controlled at the significance level that we are attempting to control.

Real Data Analysis
To illustrate the application of the proposed method to real spike train data, 15 neurons were simultaneously recorded from the cat MI shown in Figure 6 and analyzed. The experimental methodology that was implemented to collect the neural activity used for the following analysis was described in detail in Ghosh et al. [36]. Briefly, an adult cat was trained to perform a skilled reaching movement, using its preferred forelimb to retrieve food pellets placed between 2 upright Perspex barriers spaced 4 cm apart. After behavioral training was complete, PTFE coated Platinum-Iridium microwires were implanted into the cortex to a depth of about 1.5 mm into forelimb and hindlimb representations of MI (identified using intracortical microsimulation). Neural recordings were made as the animal performed the reaching task, and only neurons that significantly modulated their firing rate during task performance were isolated for analysis in this study.
Interspike interval, spike duration and spiking rate analyses were performed on neurons isolated from adjacent recording sites. This was done in order to rule out the possibility of the same neuron being counted more than once due to cross-talk between neighboring electrodes. Autocorrelogram, interspike interval and 'burst surprise' (using a surprise value of 3) analysis were performed on all neurons in order to identify any potentially bursting neurons in the data set (there were none) [37][38][39]. The data set includes 150,000 samples (3,000 samples/trial|50 trials) for each channel, and the total number of spikes for each neuron across all trials ranged from 613 to 5716. The sampling rate was 1 KHz.
Using the AIC, an optimum model for each neuron is selected to minimize the criterion. The non-overlapping spike counting window W was intuitively set to 3 ms to obtain a relatively small number of parameters while maintaining the temporal resolution.   Figure 7 shows the selected GLM spike order M i of each neuron for i~1,:::,15, and for each neuron 1 parameter (3 ms) through to 18 parameters (54 ms) were used to model its interconnection.
The GOF of the estimated model is assessed by using the Kolmogorov-Smirnov (KS) plots [40]. Prior to making inferences from an estimated statistical model, it is crucial to measure the agreement between a statistical model and the spike train data. For continuous-valued data, the GOF of the model can be quantitatively measured as standard distance such as MSE. However, this distance measure can not be applied to neural spike train data. To address this problem, we utilized the previously proposed time-rescaling theorem to transform point process measures such as neural spike train data to a continuous measure appropriate for a GOF assessment [40]. Once a CIF is estimated, rescaled times can be computed using the estimated CIF. These rescaled times will be uniformly distributed random variables on the interval (0,1 if the estimated CIF is a good approximation to the true conditional intensity of the point process. To evaluate whether the rescaled times follow the uniform distribution, we order these rescaled times from the smallest to the largest, and then plot the quantiles of the cumulative distribution function of the uniform distribution on (0,1 against the ordered rescaled times. This form of graphical representation is termed a KS plot. If the model is consistent with the data, then the points should lie on a 45-degree line. Approximate 95% confidence bounds for the degree of agreement between the model and the data may be constructed using the distribution of the KS statistic [41]. Figure 8 shows the best and the worst KS plots obtained using estimated GLMs across all the given spike train data. Most KS plots were almost within the confidence intervals, which indicates that most estimated GLMs fit the data well. The causal connectivity between the recorded neural spike train data was assessed using the proposed framework, and the results are illustrated in Figure 9. Figure 9A and B show the causal connectivity maps, Y, estimated using the proposed framework without and with the FDR correction, respectively. As illustrated in Figure 9A when the multiple comparison problem was not considered, more causal connectivity was estimated; however, there was a high probability that the false rejection of the null hypotheses of the multiple comparison caused the extra causal relationships. In the present study, a for the hypothesis testing was set to 0.05. After the FDR correction for the multiple comparison problem, the incidence of interactions between the recorded neurons was sparser, and is shown in Figure 9B. In Figure 9B, neurons 2, 3, 4, 5, 6, 9, and 10 exchanged causal interactions with a handful of other neurons including themselves, neurons 7, 12, 14, and 15 showed purely self-inhibitory interactions, and finally neurons 1, 8, 11, and 13 did not receive any influence from other neurons, nor did they show signs of selfinteraction. Interestingly, neurons 5, 6, 9, and 10 appeared to display evidence of self-excitatory interactions, which is highly unusual behavior for a neuron. Interspike interval and autocorrelogram analysis were performed on these neurons in order to exclude the possibility that these interactions were occurring due to bursting behavior [37,38]. Further analysis of these neurons revealed that they also had the four highest history orders among neurons as shown in Figure 7. Figure 9C shows the connectivity map obtained using the neural spike train data recorded during a period of postural maintenance from the same recording sites in MI following completion of a satisfactory number of task trials shown in Figure 6. The data set includes 54,000 samples (3000 sampless/ trial|18 trials), and the total number of spikes for each neuron across all trials ranged from 55 to 1030. As shown in the figure, during the state of postural maintenance, most neurons did not show any evidence of significant interactions. It could be argued that the decrease in the number of detected significant interactions that were seen during the state of postural maintenance was actually related to the decreased number of spikes that were observed during this behavioral period. In order to prove that this decrease is actually related to a physiological phenomenon rather than a decreased spike count, causality analysis was performed using the first 11 trials (of a total of 50) of the 'reaching' data set, which decreased the averaged number of spikes in that set of data (522 spikes) to a similar level as the 'postural maintenance' set (520 spikes). Figure 9D illustrates the obtained causal connectivity map, and more significant interactions were still seen between neurons during reaching movement than during postural maintenance. Note that the obtained causal connectivity maps do not necessarily represent interactions as a result of direct anatomical connection, but suggests that a functional causal connectivity exists between the recorded neurons.
The estimated GLM parameters fĉ c (0) i,i,m g Mi m~1 that correspond to the self-interactions of all neurons i for i~1,:::,15 are illustrated in Figure 10. The red, blue, and green colors represent the excitatory, inhibitory, or no self-interactions, respectively. In all cases, the first parameter is always negative due to the absolute refractory period, and the remaining parameters generally have positive values for the self-excitatory interactions and negative values for the self-inhibitory interactions, respectively. In cases where no-interactions was occurring, only one negative parameter (indicated with green asterisk) existed. Neurons showing evidence of excitatory selfinteractions have the four highest history orders, and those indicating inhibitory self-interactions have higher orders than those with no self-interactions, which have only one parameter.

Discussion
We proposed a point process framework for identifying causal relationships between simultaneously recorded multiple neural spike train data. Granger causality has proven to be an effective method to test causality between signals when using the MVAR model, but to date it has been used for continuous-valued data [7][8][9][10][11]. The method described in this study represents a novel attempt to apply Granger causality to point process data. The high level of accuracy that our method displayed when predicting the nature of the interactions occurring in the simulated data set was an encouraging indication that the proposed method is sound. Furthermore, the marked disparity in incidence of interactions during movement and non-movement periods in the experimental data is in keeping with the findings of previous studies investigating interactions in MI [36]. Thus, the outcome of both our simulated and experimental data analysis provides compelling evidence that Granger causality can be successfully applied to point process data. This is an important finding, as there are currently very few techniques that assess interactions between multiple neurons as well as providing insight regarding the causal relationships that exist between them. The ability to infer causal relationships between interacting neurons provides us with important information about networks of neurons being studied with this method. A detailed understanding of the interactions occurring in ensemble activity recorded from MI may lead to improved accuracy in algorithms used to control devices such as brain-computer interfaces and neural prosthetics [20,42,43].
Other model-based methods for assessing the directional relationships between neurons have been recently developed [21,42]. These methods infer underlying interactions between neurons based on estimated model parameters, which contain the information on the dependencies between all of the recorded neurons. Thus, functional connectivity between neurons is inferred when the estimated model parameters achieve non-zero magnitude, that is, when their confidence intervals do not cross the zeromagnitude line. However, no quantitative criteria currently exists to guide users of these methods to accept or reject detected interactions when the suspected interaction is of low magnitude, or high magnitude but with wide confidence intervals for the estimated parameters. Thus, in more difficult cases where a model-based method produces uncertain results, the acceptance of a spurious interaction, or the rejection of a legitimate one, may compromise the reliability of experimental data analysis. The proposed point process framework addresses this issue by performing statistical significance tests that investigate the causal interactions based on the likelihood ratio statistic, eliminating this uncertainty. Thus, the proposed method may be of use to researchers who are having trouble quantifying some of the connections that they are detecting when using other model-based methods.
Some of the neurons included in this analysis showed no evidence of either self-interaction, or interactions with other neurons. Although these neurons also had non-zero GLM parameters for self-interaction, as indicated with green asterisks in Figure 10, their effects were tested to be statistically insignificant compared to those caused by other neurons or background firing activity. In these cases, we must consider that the hypothesis testing to evaluate the significant causal interactions depends on the FDR value that is chosen. Decreasing the FDR value means that a statistical threshold for the significance test is more strict, and would lead to a sparser causal connectivity map. Therefore, it should be noted that the inferred causal connectivity maps Y generated by this method are not absolute, and may change depending on the user's selection of the FDR value.
The identification of excitatory self-interactions for some of the analyzed neurons was an unexpected and interesting finding. Analysis of the spiking features of these neurons verified that they were not engaged in any manner of bursting behavior that may explain the self-excitation result. Based on the high history orders that were also seen in those neurons as shown in Figure 10, we infer that the self-excitation result may be caused by 'hidden' positive feedback networks, that is, networks involving neurons that were not recorded by our microwires. To support our inference, we have performed another simulation to investigate the effect of hidden feedback networks. We identified the causal interactions among ensemble spiking activity, which was synthetically generated based on the five-neuron network of Figure 11. Compared to the nineneuron network of Figure 1A, the five-neuron network of Figure   had hidden neurons 4 and 5, which composed hidden positive feedback networks together with neuron 1. So the firing activity of neuron 1 was not only dependent on the spiking activity of observed neurons 2 and 3, but also on the spiking activity of hidden neurons 4 and 5; however, only neurons 1, 2, and 3 were observable. The parameter vector for the excitatory interaction of the hidden network was set to c (h,z) i,j, : = [0 0 1 2 2 1]. The other experimental settings were all same to the previous case. We generated 100,000 samples for each neuron, and the total number of spikes for the observed neurons 1, 2, and 3 were 4247, 2606, and 2314, respectively. Due to the hidden positive feedback, neuron 1 fired more spikes than other neurons. The model orders were selected using the AIC, and the selected orders for neurons 1, 2, and 3 were 5 (10 ms history duration), 2 (4 ms), and 2 (4 ms), respectively. Neuron 1 had a relatively longer history duration than other neurons due to the hidden feedback networks. Using the proposed method, we obtained both the Granger causality map W and the causal connectivity map Y, which is illustrated in Figure 12. The estimated causality map Y matches well the original network of Figure 11 except that neuron 1 was estimated to have a selfexcitatory interaction, which was caused by the hidden positive causal interactions with neurons 4 and 5. This hidden interaction also led to the relatively long history duration of neuron 1 compared to the other neurons. This simulation supported the idea that hidden positive feedback network leads to the relatively long history duration and can change inhibitory self-interaction to excitatory one, which we could also observe in this real data analysis case. Similarly, self-inhibitory interactions, which had a relatively long history duration as shown in Figure 10, were also identified in this study, and may be the result of hidden negative feedback networks. Self-inhibitory interactions (as they are defined using this method) may be difficult to quantify in some cases, as a neuron with a very low firing rate may produce a self-inhibitory result that is similar in appearance to that which would occur due to hidden negative networks. However, the majority of the neurons in the present study that showed the evidence of self-inhibition had quite high firing rates. Thus, the inference of hidden negative feedback networks is a plausible explanation in these cases. The proposed framework creates an unprecedented opportunity to investigate interactions from hidden neural networks that have either excitatory or inhibitory causal influences on recorded neurons. Recently a method called partial Granger causality to identify the underlying causal interactions in the presence of exogenous inputs and latent variables for the continuous-valued case has been proposed [44,45]. It would be useful to extend this work to neural spike train data in order to deal with the effects of exogenous inputs or hidden neurons beyond the investigation of the hidden feedback network.
The Matlab software and the data sets used to implement the methods presented here are available at the website (http://www. neurostat.mit.edu/gcpp).