Abstract
The functional connectivity between cortical neurons is not static and is known to exhibit contextual modulations in terms of the coupling strength. Here we hypothesized that the information flow in a cortical local circuit exhibits complex forward-and-back dynamics, and conducted Granger causality analysis between the neuronal spike trains that were simultaneously recorded from macaque inferior temporal (IT) cortex while the animals performed a visual object discrimination task. Spikes from neuron pairs with a displaced peak on the cross-correlogram (CCG) showed Granger causality in the gamma-frequency range (30–80 Hz) with the dominance in the direction consistent with the CCG peak (forward direction). Although, in a classical view, the displaced CCG peak has been interpreted as an indicative of a pauci-synaptic serial linkage, temporal dynamics of the gamma Granger causality after stimulus onset exhibited a more complex triphasic pattern, with a transient forward component followed by a slowly developing backward component and subsequent reappearance of the forward component. These triphasic dynamics of causality were not explained by the firing rate dynamics and were not observed for cell pairs that exhibited a center peak on the CCG. Furthermore, temporal dynamics of Granger causality depended on the feature configuration within the presented object. Together, these results demonstrate that the classical view of functional connectivity could be expanded to incorporate more complex forward-and-back dynamics and also imply that multistage processing in the recognition of visual objects might be implemented by multiphasic dynamics of directional information flow between single neurons in a local circuit in the IT cortex.
Introduction
Functional connectivity between single neurons has been widely investigated using cross-correlation analysis of spike trains. A center peak on the cross-correlogram (CCG) suggests common inputs or recurrent connections, whereas a displaced peak suggests pauci-synaptic serial connections (Perkel et al., 1967; Moore et al., 1970). These functional connectivities dynamically change depending on the behavioral context or presented stimulus (Engel et al., 2001; Salinas and Sejnowski, 2001; Hirabayashi and Miyashita, 2005; Kohn and Smith, 2005; Fujisawa et al., 2008).
Neuronal interactions have also been examined in the frequency domain using coherence measure (Averbeck and Lee, 2004; Fries, 2005). In particular, interactions in the gamma-frequency range have been shown to involve various cognitive processes, including attention, working memory, and object recognition (Lee, 2003; Fries et al., 2007; Zhou et al., 2008). Although coherence measure estimates the net strength of an interaction, it does not directly provide the information on the strength of the flow in each direction [but see Pesaran et al. (2008) for estimating directed interactions via spike–local field potential (LFP) coherence].
Causality analysis, including partial directed coherence, directed transfer function, and Granger causality, recently has been used to investigate the directionality of information flow for continuous neuronal signals such as LFP, EEG, or blood oxygenation level-dependent signal acquired by functional magnetic resonance imaging (Kaminski and Blinowska, 1991; Brovelli et al., 2004; Roebroeck et al., 2005; Ding et al., 2006). Although these analyses require the stationarity of the signal, which is often violated in the neuronal signals during execution of a cognitive task, several methods have been devised to cope with this issue (Ding et al., 2000; Salazar et al., 2004; Wang et al., 2008; Kayser and Logothetis, 2009). By using these methods, the dynamics of neuronal interactions between interconnected areas have been examined in several studies, the results of which suggest that information flow in each direction does not necessarily occur simultaneously but, instead, relatively serial dynamics of interdependence can be observed (Baccala and Sameshima, 1998; Kocsis and Kaminski, 2006). For example, Pesaran et al. (2008) and Gregoriou et al. (2009) demonstrated that, during a cognitive task, a top-down signal from the prefrontal cortex precedes a signal that flows in the opposite direction. Although these forward-and-back dynamics of directional interdependence have been reported in studies of macroscopic, across-areal interactions via LFP, the dynamics of directional interdependence between single neurons in a microscopic, cortical local circuit remain unknown.
In a previous study, we demonstrated that the strength of spike-timing correlation between single neurons in macaque inferior temporal (IT) cortex was modulated depending on the presented stimulus (Hirabayashi and Miyashita, 2005). Here we conducted nonparametric Granger causality analysis (Dhamala et al., 2008a,b), which is directly applicable to spike trains (Nedungadi et al., 2009), to investigate dynamic changes in the directionality of information flow between pairs of single IT neurons (Miyashita, 2004). We found triphasic temporal dynamics of directional information flow, in which the forward-and-back dynamics was followed by reappearance of the forward component. These dynamics of information flow appeared in a stimulus-dependent manner.
Materials and Methods
Behavioral tasks and recording procedures.
The experimental procedures have been described previously (Naya et al., 2003; Hirabayashi and Miyashita, 2005). In brief, two monkeys (Macaca fuscata) were trained to perform a visual object discrimination task [face-like object (FO)/nonface-like object (NFO) judgment task]. While a monkey fixed its gaze within 1.0–1.4°, a visual object (see below) was presented for 1 s (cue period). Eye position was monitored using a scleral search coil (Judge et al., 1980; Tomita et al., 1999). After a delay period (500 ms), the go signal appeared, which required the monkey to push the right or left button within 1 s depending on whether the presented stimulus was FO or NFO. All the analyzed neuronal signals in the present study were obtained during this task. Monkeys performed this task with 99.7 ± 0.4% (mean ± SEM) correct responses during recording sessions. Visual stimuli were presented on the monitor (refresh rate, 60 Hz) by using a custom-made software written in C. In a typical recording session, monkeys performed over 1000 trials per day. The presented stimulus was a whole object that was composed of four individual parts arranged in facial (FO) or random (NFO) configurations within a radius of ∼3° against a surrounding contour (7.8° high and 6.1° wide). We prepared 120 facial-feature-like parts (40 eye-like parts, 40 nose-like parts, and 40 mouth-like parts) to construct 403 = 64,000 whole objects. Before the training of the above discrimination task, monkeys were well familiarized with these 120 parts using a delayed-matching-to-sample task.
Multiple single-unit activities (SUAs) were recorded from the IT cortex in three hemispheres of the two monkeys using a tetrode (Thomas Recording) (Aronov et al., 2003; Soteropoulos and Baker, 2006). Recorded signals were amplified (JH-440J and AB-651J; Nihon Kohden) and bandpass filtered (500 Hz to 5 kHz) (Tamura et al., 2004; Fujisawa et al., 2008). During a recording session, a pair of single units was sorted using the window discrimination technique to determine the optimal whole objects for the recorded cell pair, and neuronal signals were also stored and digitized off-line at 25 kHz to sort into multiple single units by waveform analysis (see below). Only these off-line-sorted single units were further analyzed in the present study.
The parts comprising the optimal whole objects [optimal FO (pFO) and optimal NFO (pNFO)] were determined by the minimax algorithm (Baky et al., 1981; Hirabayashi and Miyashita, 2005) for each on-line-sorted neuron pair so that they elicited as high a response as possible from both neurons. For a given pair of neurons, the same set of parts was used for its pFO and the corresponding pNFO. The spatial arrangement of the parts in the pNFO was changed in a different recording session, and the same arrangement was never used again. Optimal and nonoptimal FOs/NFOs were presented to the animal in a pseudorandom order. Neuronal correlation was not examined during recording sessions, ensuring that the recorded cell pairs or their optimal stimuli could not be selected so that the cell pairs revealed spike correlation in response to their optimal whole objects. All animal experiments were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and the regulations of the University of Tokyo School of Medicine.
Off-line spike sorting.
We used the same database as our previous study (Hirabayashi and Miyashita, 2005). In the present study, the methods for spike sorting were modified as follows. We conducted off-line spike sorting using Klustakwik (Harris et al., 2000), an open source for semiautomatic spike sorting, and by following manual adjustment of the clusters (Barthó et al., 2004; Fujisawa et al., 2008). We also conducted the spike sorting by means of ICSort, the automatic spike-sorting algorithm using the independent component analysis, enabling us to sort the temporally overlapped spike waveforms into multiple single units (Takahashi et al., 2003; Sakurai et al., 2004; Sakurai and Takahashi, 2006). If the results obtained by ICSort were qualitatively similar to those obtained by Klustakwik, except for overlapping spikes, we adopted the results obtained by ICSort. In case that single units were sorted by Klustakwik, we obtained CCGs with a center dip at the bins within a 1 ms lag because of the inability to sort the temporally overlapped spike waveforms into their constituting multiple single units (Barthó et al., 2004; Tamura et al., 2004; Fujisawa et al., 2008). In such cases, we did not analyze the peak structure of the CCG at ±1 ms of the zero lag point (Barthó et al., 2004; Tamura et al., 2004; Fujisawa et al., 2008). This might have resulted in exclusion of the narrow center peak that extended within ±1 ms of the zero lag point. As for displaced peaks, however, the peaks should be detected, in most cases, with the lag time larger than 1 ms because of conduction delay between neurons, and thus we considered that displaced peaks were not excluded using this method. Monosynaptic or polysynaptic serial connections have also been examined in previous studies by detection of displaced CCG peak without bins within a 1 ms lag (Menz and Freeman, 2003, 2004; Tamura et al., 2004; Maurer et al., 2006; Fujisawa et al., 2008).
Cross-correlation analysis.
The neuronal data were analyzed only for the correct trials. We defined a cell as responsive to a whole object if the firing rate during the cue period (assessed in a 500 ms period from 80 ms after cue onset) was significantly (p < 0.05, paired t test) different from that in the corresponding period just before cue onset. We conducted cross-correlation analysis for a cell pair only if both constituent cells showed significant responses to both their pFO and pNFO that were determined during the recording session before the cross-correlation analysis.
We constructed CCGs for lag times within 100 ms with a 1 ms resolution using spikes recorded during a 1 s period beginning 80 ms after cue onset. A CCG for the pFO or pNFO was constructed from spikes recorded in at least 55 trials for each stimulus and was accepted only if the available spikes exceeded 1900 (mean ± SEM; 5022 ± 310 for all of the analyzed cell pairs, at least 700 per cell). A shift predictor, calculated using one-trial-shifted spike trains (Constantinidis et al., 2001, 2002; Tamura et al., 2004; Kohn and Smith, 2005), was smoothed (five-bin boxcar averaging) (Nowak et al., 1995) and subtracted from the raw CCG to remove the stimulus-locked component (Perkel et al., 1967; Nowak et al., 1995; Kohn and Smith, 2005), yielding the shift predictor-subtracted CCG (SSCC). The peak height of an SSCC was normalized by the SD of the shift predictor to calculate a z score (Constantinidis et al., 2001, 2002) and was identified as significant when the z score exceeded the threshold corresponding to p = 0.05 (one-sided, z > 2.81, detected within the lag time of 10 ms, corrected for multiple comparisons) (Hirabayashi and Miyashita, 2005). SSCCs were constructed from spikes recorded in 145 ± 6 (mean ± SEM) trials for each stimulus. Significant SSCC peaks were categorized into center or displaced peaks as follows. First, we calculated the asymmetry index (AI) of the peak structure of each SSCC (Alonso and Martinez, 1998; Menz and Freeman, 2003, 2004) as follows:
where R and L indicate the summed bin counts on the right and left side of the SSCC within 10 ms lags, respectively. Bins were included in the calculation only if the z score exceeded 2. A given significant SSCC was defined to possess a displaced peak if the absolute value of the AI was ≥0.5, whereas it was defined to possess a center peak if the absolute value of the AI was <0.5. A displaced peak on the SSCC suggests pauci-synaptic serial connections, whereas a center peak suggests common inputs or recurrent connections (Perkel et al., 1967; Moore et al., 1970). Although the above definition is somewhat arbitrary as suggested in previous studies (Menz and Freeman, 2003, 2004), still it reflects the differences in the underlying anatomy of the observed functional connectivity and thus is effective for separately investigating their functional properties, as is the case in the present study and in previous studies (Alonso and Martinez, 1998; Menz and Freeman, 2003, 2004) (see Discussion for details).
Spectral analysis of spike trains.
Spectral estimates of spike trains (auto- or cross-spectrum and coherence) were calculated using the multitaper method (Jarvis and Mitra, 2001; Pesaran et al., 2002, 2008; Saalmann et al., 2007; Fries et al., 2008; Gregoriou et al., 2009; Nedungadi et al., 2009) using Chronux, an open-source MATLAB software package (http://www.chronux.org) (DeCoteau et al., 2007; Zhou et al., 2008). In this method, spectral estimates of spike trains were determined as the average of the results that were computed after multiplying several different taper functions to the original signal. These processes yield less-biased and noise-reduced spectral estimates compared with conventional methods (Mitra and Pesaran, 1999; Jarvis and Mitra, 2001). For time-averaged analysis in which spectral estimates were calculated for 1-s-long spike trains, we applied nine orthogonal Slepian tapers to calculate the spectral estimates of spike trains that were downsampled to 1 kHz, providing a resultant frequency resolution of ±5 Hz (Pesaran et al., 2002). In time-resolved analysis, spectral estimates of spike trains were calculated by using five-taper functions for a ±125 ms window that was slid in 25 ms steps (frequency resolution, ±12 Hz) (Zhou et al., 2008).
Nonparametric Granger causality analysis.
The Granger causality estimates whether the value of a given time series is better predicted by incorporating the past knowledge of another time series compared with the case that the prediction was conducted by using only the past knowledge of the first time series itself. Ordinary parametric Granger causality (Brovelli et al., 2004; Ding et al., 2006; Gregoriou et al., 2009; Kayser and Logothetis, 2009) and some other analyses for estimating directional interactions (Kaminski and Blinowska, 1991; Baccala and Sameshima, 1998) require multivariate autoregressive (MVAR) modeling of the signal. To apply these methods to spike trains, they have to be converted into continuous data, which might disrupt the properties of spike trains as point processes (Truccolo et al., 2005; Dhamala et al., 2008a,b; Nedungadi et al., 2009). In the present study, Granger causality analysis was conducted using a nonparametric method, which bypasses MVAR modeling of the data (Dhamala et al., 2008a,b), and is thus directly applicable to spike trains (Nedungadi et al., 2009).
We conducted nonparametric Granger causality analysis according to the method described in previous studies (Dhamala et al., 2008a,b; Nedungadi et al., 2009). In brief, spectral estimates of spike trains were first calculated by means of the multitaper method, as explained above. To calculate the Granger causality, transfer function and noise covariance matrix were estimated from the spectral estimate S(f) using spectral matrix factorization (Wilson, 1972), which decomposes S(f) into a unique corresponding transfer function H(f) and noise covariance matrix V (Dhamala et al., 2008a,b; Nedungadi et al., 2009). According to Geweke's formulation (Geweke, 1982, 1984), the auto-spectrum of a given signal at a given frequency can be decomposed into an intrinsic component and the other component that is predicted from another signal. The ratio of the auto-spectrum to the intrinsic component thus provides the amplitude of the causal influence. Here, the Granger causality from spike train Nj(t) to Ni(t) at frequency f can be represented as follows:
The Granger causality in the opposite direction can be estimated similarly. Spectral factorization and estimation of the Granger causality were conducted by the MATLAB scripts that were kindly provided by Prof. G. Rangarajan (Department of Mathematics, Indian Institute of Science, Bangalore, India) (Nedungadi et al., 2009). For a given cell pair, coherence and Granger causality values in Figures 1⇓⇓⇓⇓–6 were calculated using spike trains for the pFO when the SSCC exhibited a significant displaced peak for both the pFO and pNFO. When only the pNFO elicited a significant displaced peak, spike trains for the pNFO were used for these coherence and Granger causality analyses. First, we examined the time-averaged Granger causality in the cue period by using spike trains in a 1 s period from 80 ms after cue onset. Next, we conducted time-resolved analysis of the Granger causality (Ding et al., 2000; Zhang et al., 2008; Gregoriou et al., 2009) with a ±125 ms window that was slid in 25 ms steps. Note that in this time-resolved analysis, stationarity of the signal is expected to be improved because of segmentation of the signal into a short window when compared with the time-averaged analysis described above (Ding et al., 2000; Gregoriou et al., 2009) (see Discussion for details).
Statistical test for significance of coherence and Granger causality.
We adopted a nonparametric permutation test to examine statistical significance of coherence for each cell pair (Maris et al., 2007; Fries et al., 2008). The null hypothesis was that the obtained value of coherence was not significantly different from those calculated with trial-shuffled spike trains. For a given pair of neurons, coherence was calculated for trial-shuffled spike trains, and the maximum value of coherence across all the frequencies in the gamma range (30–100 Hz) was determined. These procedures were repeated 1000 times to obtain a null distribution of the maximum gamma coherence. The experimental value of gamma coherence was considered to be significant if it exceeded the 50th value from the maximum of the aforementioned null distribution, which was equivalent to a p value of 0.05 for a one-sided test. The same procedure was also applied for the Granger causality in each direction. To test the significance of directionality, the threshold was set at the 25th value from the maximum of the null distribution, which corresponded to a p value of 0.05 for a two-sided test. This method provided a statistical evaluation corrected for multiple comparisons across different frequencies in the gamma range, because all the values constituting the null distribution were determined as the maximum across different frequencies (Maris et al., 2007).
The randomization procedure was also applied for the population analysis of time-resolved data to correct for multiple comparisons across different time points (Womelsdorf et al., 2007; Fries et al., 2008). In that case, trial-shifted coherence or Granger causality was subtracted from the original data for each cell pair before testing the statistical significance, and thus the null hypothesis was that the population average of the resultant value was not significantly larger than zero. To estimate the null distribution, the veridical data set was first randomly divided into two groups, and the mean time course was calculated for each group. After subtraction of one time course from the other, the maximum value was determined across all the time points. These processes were repeated 1000 times to yield the null distribution, to which the population average time course obtained from the experiment was compared. The threshold was set at the value corresponding to a p value of 0.05 for a one-sided test. For testing the significance of directionality, the threshold was set at the value corresponding to a p value of 0.05 for a two-sided test.
Parametric Granger causality analysis.
To apply the MVAR modeling to spike trains during the cue period (a 1 s period from 80 ms after cue onset), we first converted spike trains into continuous data (Jarvis and Mitra, 2001; Nedungadi et al., 2009). In brief, cumulative spike count of a given neuron, N(t), at time t was converted into a zero-mean continuous signal, N̅(t), by subtracting the expected spike count λt:
where λ is the mean firing rate of the neuron. After this conversion of the spike trains, signals in each trial were sampled at 200 Hz. From this point, we conducted parametric Granger causality analysis by using the open-source MATLAB package (BSMART) (Cui et al., 2008). In brief, the converted neuronal signals were first normalized by subtracting the point-by-point ensemble mean across trials and dividing by the point-by-point SD across trials. These processes eliminated nonstationarity of the signal that stems from the temporal dynamics of the mean and SD. The above data were then fitted to the bivariate autoregressive model, in which a pair of above signals X(t) as a function of time t was written as follows:
where E(t) is the two-dimensional noise vector as a function of time t, Am is the 2 × 2 coefficient matrix, and p represents the model order. To determine the model order of MVAR modeling, Akaike's information criteria (AIC) was calculated for each pair of neurons. As a population, the minimum AIC value for cell pairs with a displaced CCG peak was 6.7 ± 0.49 (mean ± SEM; n = 21), and thus we adopted a model order of 7 to conduct the MVAR modeling of our data. The transfer function H(f) was computed from the coefficient matrix Am as follows:
The estimated spectral matrix was given by the following:
where * denotes matrix transposition and complex conjugation and V is the covariance matrix of the noise vector. Once the spectral estimate was computed, Granger causality in the frequency domain was estimated by using Geweke's decomposition (Geweke, 1982, 1984).
Analysis of MS-triggered neuronal responses.
Horizontal and vertical eye traces were sampled at 10 kHz using Recorder (Plexon) and were downsampled to 1 kHz. The velocity of eye movement was calculated as a square root of the sum of squares of deviations in vertical and horizontal directions compared with those at the previous time point. Detection of microsaccades (MSs) was conducted basically in accord with the method described in a previous study (Martinez-Conde et al., 2000). MSs were detected in the cue period during which the monkeys fixated within a window of ±1.4°. Eye movement was defined as MS if the velocity exceeded 5°/s and if the displacement exceeded 3 arcmin. The effects of MSs on the neuronal responses were examined primarily based on the analyses described previously (Bosman et al., 2009). First, the MS-triggered average of the firing rate for a given neuron was calculated as follows: for each MS event in the cue period during which the pFO or pNFO was presented, the firing rate in a time window of ±500 ms around the MS onset was calculated with a resolution of 10 ms. This histogram was calculated for all MS events in the cue period, and the average histogram was labeled as the raw histogram of the MS-triggered average of the firing rate. To control for the effect of the stimulus-locked response, we calculated the same histogram after trial shifting of spike trains against MS events. This histogram was subtracted from the raw histogram to obtain the shift-subtracted MS-triggered average of the firing rate. To evaluate the practical impact of MSs on the actual neuronal responses during a trial, we compiled the time course of MS-induced modulation of the firing rates by convolving the timing of MS onset (10 ms resolution) with the shift-subtracted MS-triggered average of the firing rate. The convolution was conducted for all trials, and the resultant time course was normalized by the number of trials. Similar procedures were repeated to calculate the effect of MSs on the coherence and Granger causality, except that these were calculated for pairs of neurons using a ±125 ms sliding window that was slid in 25 ms steps. We adopted a nonparametric permutation test (Maris et al., 2007; Womelsdorf et al., 2007; Fries et al., 2008) to examine the statistical significance of the effect of MSs on the neuronal responses. The statistical significance of the shift-subtracted MS-triggered average and the MS-induced modulation of coherence and Granger causality were assessed for the gamma component, as in the original methods for examining the statistical significance of stimulus-elicited coherence and Granger causality (see Fig. 5A,B).
All of the statistical tests in the present study were two-sided unless stated otherwise.
Results
Gamma Granger causality between spike trains of single IT neurons
Simultaneous recordings were conducted from multiple single units in macaque IT cortex, while the animals performed a visual discrimination task. We recorded the activity of 92 responsive cells with a sufficient number of spikes in the cue period for cross-correlation analysis (see Materials and Methods; 28 and 64 cells from monkeys 1 and 2, respectively). These cells constituted 61 cell pairs in total (20 and 41 cell pairs for monkeys 1 and 2, respectively). Of these, 21 cell pairs (12 and 9 cell pairs for monkeys 1 and 2, respectively) revealed a significant displaced peak on the CCG in response to their optimal stimulus [see Materials and Methods; peak z score, 5.0 ± 0.88 (mean ± SEM); peak lag, 3.6 ± 0.47 ms; AI, 0.90 ± 0.04; half-width at half-maximum, 1.8 ± 0.33 ms]. In the present study, we analyzed the spectral properties of spike correlation for these cell pairs exhibiting a displaced CCG peak.
The representative data of the spectral analysis for a pair of IT neurons is depicted in Figure 1. Both neurons of this pair exhibited significant spiking responses to their optimal stimulus (Fig. 1A). For this cell pair, we calculated the CCG between the spike trains during the cue period and found that this cell pair exhibited a significant displaced peak on the left side of the CCG (peak z score, 3.9; AI, −1.0; peak lag, −1.5 ms; half-width at half-maximum, 3.5 ms) (Fig. 1B), suggesting pauci-synaptic connections from unit 1 to unit 2. Next, we applied multitaper spectral analysis to this cell pair for investigating the spectral characteristics of the functional connectivity (Jarvis and Mitra, 2001; Pesaran et al., 2002, 2008; Saalmann et al., 2007; Fries et al., 2008; Zhou et al., 2008; Gregoriou et al., 2009). Spike trains of this cell pair revealed prominent gamma coherence during the cue period (Fig. 1C, orange trace). Trial shifting (Saalmann et al., 2007; Zhou et al., 2008; Gregoriou et al., 2009) essentially diminished this gamma coherence (Fig. 1C, khaki trace), suggesting that the observed coherence was not attributable to the stimulus-locked dynamics of spike firings. Statistical significance of the coherence was evaluated using a null distribution of the data calculated from trial-shuffled spike trains (Maris et al., 2007), and the peak value of the gamma coherence of this cell pair was found to be significant (peak, 41 Hz; p < 0.001, permutation test, corrected for multiple comparisons across frequencies). The results of nonparametric Granger causality analysis for the above cell pair are depicted in Figure 1D. In both directions, Granger causality in the gamma range revealed prominent peaks at the frequency similar to that of coherence, but the peak value in the direction from unit 1 to unit 2 was greater than that in the opposite direction (Fig. 1D, red and blue traces), consistent with the results of cross-correlation analysis. Trial shifting (Brovelli et al., 2004; Gregoriou et al., 2009; Nedungadi et al., 2009) almost completely attenuated the causality in both directions (Fig. 1D, magenta and cyan traces). Statistical significance of the gamma Granger causality was evaluated using trial-shuffled spike trains as in the case of coherence (Maris et al., 2007), and the peak value of the Granger causality was significant (unit 1 to unit 2: peak, 38 Hz; p < 0.001, permutation test; unit 2 to unit 1: peak, 42 Hz; p < 0.002). Directionality of the Granger causality was defined as the difference in the causality value between the two directions (unit 1 to unit 2′ minus unit 2 to unit 1′), and the directionality for this cell pair was found to be significant (black trace; peak, 38 Hz; p < 0.001, permutation test) (Fig. 1E).
Functional connectivity between a pair of IT neurons. A, Responses of a pair of single IT neurons to their optimal stimulus. The horizontal bar below the histogram indicates 1 s duration of stimulus presentation (cue period). B, CCG between the spike trains during the cue period. Left, Raw CCG (black) and shift predictor (gray). Abscissa, Spike time of unit 1 relative to that of unit 2. Bin width, 1 ms. Right, SSCC of the cell pair. Horizontal gray lines, Confidence limit (p = 0.05, corrected for multiple comparisons). C, Coherence between raw spike trains (orange) and between trial-shifted spike trains (kahki). Horizontal orange line, Confidence limit (p = 0.05, corrected for multiple comparisons). D, Granger causality between raw spike trains (red, unit 1 to unit 2; blue, unit 2 to unit 1) and between trial-shifted spike trains (magenta, unit 1 to unit 2; cyan, unit 2 to unit 1). Red and blue horizontal lines, Confidence limit for unit 1 to unit 2 (red) and unit 2 to unit 1 (blue). E, Directionality of the Granger causality (unit 1 to unit 2 minus unit 2 to unit 1) for raw (black) and trial-shifted (gray) spike trains. Horizontal black line, Confidence limit.
In the population of cell pairs with the displaced CCG peak (n = 21), gamma coherence in the cue period was greater than the control that was calculated using trial-shifted spike trains (peak, 51 Hz) (Fig. 2A, orange). In the precue baseline period (1000–0 ms before cue onset), coherence was much weaker than that in the cue period (Fig. 2D, orange). Two-way ANOVA for gamma coherence (30–80 Hz) showed that gamma coherence in the cue period was significantly greater than the trial-shifted control or than the gamma coherence in the baseline period [p < 0.001, two-way ANOVA followed by least significant difference (LSD); n = 21] (Fig. 2G). Significant gamma coherence in the cue period was consistently observed in both monkeys (monkey 1: p < 0.001, paired t test, n = 12; monkey 2: p < 0.004, n = 9).
Population data of coherence and Granger causality for cell pairs with a displaced CCG peak. A–C, Coherence and Granger causality in the cue period. A, Coherence for raw (orange) and trial-shifted (kahki) spike trains (n = 21). Thick and thin traces depict the mean and mean ± SEM, respectively. B, Granger causality for raw (red and blue) and trial-shifted (magenta and cyan) spike trains. The direction of causality was defined as forward and backward according to the CCG peak (see Results). C, Directionality of Granger causality for raw (black) and trial-shifted (gray) spike trains. The directionality was defined as the signed difference between the forward and backward components of Granger causality. D–F, Same analyses as in A–C, but for spike trains during the baseline period (1000–0 ms before cue onset). G, H, Average gamma (30–80 Hz) components of coherence (G) and Granger causality (H) for the population. Error bars represent SEM. Asterisks indicate statistical significance assessed by two-way ANOVA followed by LSD.
Population coherence spectra revealed a broader peak compared with the representative data in Figure 1. Even in the case that the coherence of each cell pair had a clear frequency structure stemming from an oscillatory process, the population data of the coherence spectrum might become broader than that from a single-cell pair, if the peak frequency of each coherence varied pair by pair. To assess this possibility, we examined the peak frequency and frequency width (FWHM) of each coherence (supplemental Fig. 1, available at www.jneurosci.org as supplemental material). The median value of the FWHM was 23.7 Hz, and more than half (57%) of the cell pairs revealed coherence with the FWHM that was smaller than 30 Hz (supplemental Fig. 1A, available at www.jneurosci.org as supplemental material). The peak frequency of coherence varied cell pair by cell pair, and the median value was 54.2 Hz (supplemental Fig. 1B, available at www.jneurosci.org as supplemental material). These results indicate that, as described above, individual spike–spike coherences in the present study revealed distinct peaks with various peak frequencies in the gamma range, suggesting the involvement of underlying oscillatory processes.
Next, we examined the phase spectrum of coherence for the frequency range in which significant coherence was observed with the permutation test. Figure 3, A and B, shows the superposed phase spectra of coherence for individual cell pairs that were calculated using raw and trial-shifted spike trains, respectively. In population, the phase of coherence for raw spike trains was significantly more stable (i.e., with smaller variance) across frequencies than that which was calculated using trial-shifted spike trains (p < 0.001, paired t test). These results indicate that pairs of single IT neurons with a displaced CCG peak exhibited reliable gamma coherence with phase stability across different frequencies.
Phase of coherence and its relationship with the directionality of Granger causality. A, B, Superposed phase spectra of coherence for individual cell pairs with a displaced CCG peak. Phase spectra of coherence were calculated using raw (A) and trial-shifted (B) spike trains. Different colors correspond to different cell pairs. Phase of coherence was calculated for frequencies in which the amplitude of coherence was significant for raw spike trains and the mean phase within the calculated frequencies was subtracted. Note that 0 Hz depicts the lowest frequency at which coherence of a given cell pair was significant. C, D, Relationship between the phase of coherence and the directionality index of Granger causality. Granger causality was calculated using raw (C) and trial-shifted (D) spike trains. E, Relationship between the phase of coherence and the directionality index of Granger causality calculated for individual frequencies. Horizontal line, Confidence limit (p = 0.05, corrected for multiple comparisons).
For the population analysis of the Granger causality, directions of the influence for each cell pair were assigned according to the results of cross-correlation analysis: when a given CCG revealed a significant displaced peak indicating the functional connectivity from unit 1 to unit 2, the direction from unit 1 to unit 2 was defined as “forward,” whereas the opposite was defined as “backward.” As a population, the gamma Granger causality in the forward direction during the cue period was greater than the trial-shifted control or the corresponding value in the backward direction (peak, 47 Hz) (Fig. 2B, red). Directionality of the Granger causality revealed a prominent peak in the gamma range (peak, 48 Hz) (Fig. 2C, black). During the cue period, the gamma Granger causality (30–80 Hz) in either direction was significantly larger than the trial-shifted control (forward, p < 0.01; backward, p < 0.04) (Fig. 2H, left), and the value in the forward direction was significantly larger than that in the backward direction (p < 0.02) (Fig. 2H, left). Forward-dominant causality in the gamma range was consistently observed in both monkeys [monkey 1: forward/backward, 3.01 ± 0.87 (mean ± SEM; p < 0.03, one-sided paired t test; n = 12); monkey 2: 2.94 ± 0.89 (p < 0.04; n = 9)]. In the baseline period, however, the gamma Granger causality in either direction was not significantly different from the trial-shifted control, and the values in each direction were not significantly different from each other (p > 0.07) (Fig. 2E,F,H right). Together, these results indicate that pairs of single neurons in the IT cortex with the displaced CCG peak exhibit Granger causality in the gamma range in response to the presented stimulus, with a directionality in accordance with the CCG peak.
In the present study, Granger causality was calculated using a nonparametric method, which bypasses MVAR modeling of the data (Dhamala et al., 2008a,b), and it enables direct application to spike trains (Nedungadi et al., 2009). To confirm that the present results do not depend on the adopted method for estimating the Granger causality, we also conducted parametric Granger causality analysis (see Materials and Methods). For the spike trains that exhibited a displaced CCG peak, the parametric Granger causality analysis in the cue period showed qualitatively similar results to those obtained using the nonparametric method: forward causality showed a prominent peak in the gamma range (49 Hz), and the gamma component of forward causality (30–80 Hz) was significantly greater than that in the backward direction (p < 0.05, two-way ANOVA followed by LSD; n = 21) or that calculated with the trial-shifted spike trains (p < 0.007). These data confirm that our results were not attributable to a specific method for estimating the Granger causality.
Next, we examined whether the strength of directionality of the Granger causality corresponded with that of the CCG. The directionality index of the Granger causality was defined as the normalized difference between the forward and backward components. The directionality index of the gamma Granger causality (30–80 Hz) was significantly correlated with the AI of the corresponding CCG (r = 0.652; p < 0.001) (supplemental Fig. 2A, available at www.jneurosci.org as supplemental material). This relationship was not observed when the trials were shifted for calculating the Granger causality (r = −0.201; p > 0.3) (supplemental Fig. 2B, available at www.jneurosci.org as supplemental material). We also examined whether the directionality of the Granger causality is correlated with the phase of coherence, which reflects the systematic lag time between the correlated signals (Kocsis and Vertes, 1997; Schack et al., 1999; Bressler and Kelso, 2001; Weiss and Mueller, 2003; DeCoteau et al., 2007; Saalmann et al., 2007; Gregoriou et al., 2009). The directionality index of the Granger causality was also found to be significantly correlated with the phase of coherence (r = 0.577; p < 0.002) (Fig. 3C). Trial shifting for calculating the Granger causality, again, essentially eliminated this correlation (r = −0.215; p > 0.2) (Fig. 3D). These results suggest that for a given cell pair, the directionality of the gamma Granger causality is closely related to both the asymmetry of the CCG peak and the phase of coherence.
To further examine the above relationship in more detail, we calculated the relationships between the phase of coherence and the directionality of Granger causality for individual frequencies. We conducted the calculation with steps of 10 Hz to obtain the results practically for individual frequencies considering the analytical resolution (see Materials and Methods). The obtained results are displayed in supplemental Figure 3 (available at www.jneurosci.org as supplemental material), and are summarized as a function of the frequency in Figure 3E. Statistically significant correlations between the phase of coherence and the directionality index of Granger causality were found for the frequency ranges of 40–50 and 50–60 Hz (p < 0.05 and p < 0.02, respectively, corrected for multiple comparisons). These results indicate that a significant positive correlation between the phase of coherence and the directionality of Granger causality was observed even when the relationships were examined for individual frequencies in the gamma range. Note that coherence and Granger causality were most prominent in these frequency ranges (Fig. 2A,B).
Triphasic rapid dynamics of directional interdependence between single neurons
The above Granger causality analysis evaluated a temporal average value of directional influences in the cue period. We next examined the dynamics of the interdependence using the sliding window technique (Ding et al., 2000; Zhang et al., 2008; Gregoriou et al., 2009). We calculated the Granger causality in a ±125 ms window with 25 ms steps and subtracted those calculations with trial-shifted spike trains. The results of a representative cell pair are depicted in Figure 4A–E. This cell pair exhibited gamma coherence throughout the cue period (peak frequency, 46 Hz; peak latency, 225 ms from cue onset) (Fig. 4A). Granger causality of this cell pair also revealed a transient increase with the dominance in the forward direction (peak frequency, 50 Hz; peak latency, 250 ms) (Fig. 4B). In a later epoch, a transient forward causality reappeared (peak frequency, 58 Hz; peak latency, 550 ms). Directionality analysis showed that forward-dominant transient causality was observed twice in the cue period (Fig. 4D). Firing rate dynamics of this cell pair did not reveal the corresponding double-peak structure (Fig. 4E), and the early peak of both coherence and Granger causality appeared before the peak of spike firing.
Temporal dynamics of coherence and Granger causality. A–E, Data from a pair of neurons. A, Temporal dynamics of coherence calculated by sliding a ±125 ms window in 25 ms steps. Time 0 depicts the window in which spikes from −125 to 125 ms from cue onset were included. Data obtained using trial-shifted spike trains were subtracted. B, C, Temporal dynamics of Granger causality in the forward (B) and backward (C) directions. D, Temporal dynamics of directionality. E, Firing rates of the two neurons. Note that the peak firing rates were different between the units (see the left and right of the color bar for units 1 and 2, respectively). F–J, Same analyses as in A–E, but for the population of cell pairs with a displaced CCG peak (n = 21). Pre-units and post-units in J were defined from the CCG peak (see Results).
The population data of this time-resolved spectral analysis for cell pairs with a displaced CCG peak (n = 21) are depicted in Figure 4F–I. Population gamma coherence revealed a sustained increase in the cue period (Fig. 4F). Population Granger causality revealed a dominance in the forward direction (Fig. 4G,H). Compared with the coherence, forward causality exhibited more transient dynamics. The early component was transiently observed in a relatively high-frequency range (peak frequency, 74 Hz; peak latency, 100 ms) and was rapidly diminished despite continuous increase in the coherence. The forward causality transiently reappeared in a later epoch of the cue period (peak frequency, 58 Hz; peak latency, 600 ms). In contrast to the forward causality, the backward causality exhibited a small and slow increase when the early component of the forward causality was attenuated, and it occurred only in a relatively low-frequency range (Fig. 4H). Directionality analysis showed two separate peaks of forward dominance [i.e., a strong early peak (peak frequency, 78 Hz; peak latency, 125 ms) and a relatively weak late peak (peak frequency, 58 Hz; peak latency, 600 ms)] (Fig. 4I). The population firing rate was examined by defining the pre-units and post-units from the direction of functional connectivity inferred from the CCG. Firing rate dynamics were very similar between the pre-units and post-units (Fig. 4J), and their peak latencies (250 ms for both pre-units and post-units) were longer than that of the early peak of the forward causality. In contrast to the dynamics of forward Granger causality, firing rates of single neurons did not reveal a substantial rebound after the first peak. These results indicate that neuronal interactions exhibit dynamic changes not only in terms of the strength (coherence) but also in terms of the directionality (Granger causality) and that the rapid dynamics of directed influences is not merely a reflection of firing rate dynamics of individual neurons.
To quantify these rapid changes in the neuronal interaction, we calculated the dynamics of the gamma component for the coherence and Granger causality. Gamma coherence (30–80 Hz) revealed a significant sustained increase throughout the cue period (peak, 275 ms from cue onset; p < 0.001, permutation test, corrected for multiple comparisons across time bins) (Womelsdorf et al., 2007; Fries et al., 2008) (Fig. 5A). Granger causality exhibited more dynamic changes: gamma Granger causality in the forward direction (30–80 Hz) revealed two significant detached peaks (150 and 550 ms; p < 0.001 and p < 0.03, respectively, permutation test as in the case of coherence) (Fig. 5B, red trace). To compare the dynamics of coherence and forward causality, we calculated the FWHM of these time courses. FWHM of the forward causality dynamics was significantly shorter than that of the coherence dynamics (193 ± 24 vs 402 ± 54 ms, mean ± SEM; p < 0.001, paired t test). Gamma Granger causality in the backward direction (30–65 Hz) showed a single, slowly developing significant peak when the early component of the forward causality was attenuated (peak: 400 ms, p < 0.004) (Fig. 5B, blue trace). We confirmed that even if the same analysis was conducted by only using cell pairs that revealed the highest AI in their CCG (AI, 1.0; n = 15), the slowly developing backward causality was still significant (peak: 400 ms, p < 0.009), suggesting that this slowly developing backward causality is not attributable to cell pairs whose CCG peak structure is relatively close to that of the center peak. The directionality dynamics of gamma Granger causality (30–80 Hz) revealed two separate peaks that corresponded with the dynamics of the forward component, although only the first peak was significant (150 ms; p < 0.02) (Fig. 5C). Next, we compared the latency of forward and backward causality. Half-maximum latency of the backward causality was significantly longer than that of the forward causality (360 ± 50 vs 239 ± 45 ms, mean ± SEM; p < 0.04, paired t test). This result suggests that the spiking activity in a local circuit was transferred in the forward direction first and then retransferred in the backward direction with a substantial delay. There were no systematic differences either in the firing rates between the pre-units and post-units throughout the cue period or in the response latencies between these groups of single neurons (Fig. 5D). Furthermore, both the pre-units and post-units showed only a single peak in their firing rate dynamics. These results suggest that the observed causality dynamics and its forward dominance was not a simple reflection of firing rate dynamics of either pre-units or post-units.
Statistical analysis of coherence and Granger causality dynamics in the gamma range. A, Time course of gamma coherence (n = 21). Thick and thin traces are mean and mean ± SEM, respectively. Horizontal orange line, Confidence limit (p = 0.05, corrected for multiple comparisons). B, Time course of gamma Granger causality. Horizontal red and blue lines, Confidence limit for the forward and backward directions, respectively. C, Time course of directionality. Horizontal black line, Confidence limit. D, Time course of firing rates for the pre-units and post-units.
We derived three significant components of the Granger causality (i.e., early forward, middle backward, and late forward components) and examined the peak frequency of each component. The peak frequency of the early forward component was significantly higher than that of the middle backward component (58.5 ± 3.8 vs 47.8 ± 3.2 Hz, mean ± SEM; p < 0.006, paired t test, corrected for multiple comparisons) (Fig. 6). The peak frequency of the late forward causality (52.2 ± 3.5 Hz) was not significantly different from either of the other two components (p > 0.2).
Peak frequency of Granger causality. Three significant components (early forward, middle backward, and late forward) were derived from the time course of gamma Granger causality. Each histogram represents the population average (n = 21) of the peak frequency for each component of Granger causality. Error bars represent SEM. The asterisk indicates statistical significance assessed by the paired t test, corrected for multiple comparisons.
In the above analyses of Granger causality dynamics for displaced CCG peak pairs, we defined the pre-unit and post-unit for each cell pair on the basis of the direction of functional connectivity inferred from the CCG. This raises the concern that the existence of backward causality might be a reflection of erroneous assignment of the direction. When we evaluate the causality of a single pair of neurons, exchanging the order of the spike trains (unit 1 or 2) leads to just reversing the triphasic pattern. However, when we conducted a population analysis of the Granger causality, it was necessary to define the order of the spike trains on the basis of a common rule to compile the data on individual cell pairs into population data. In the first analysis, we used the CCG peak to define the order of the spike trains for individual cell pairs. We could also define the order on the basis of the Granger causality itself. This other method of defining could potentially result in a different order of spike trains for the same cell pair. Especially, the backward causality, which is generally weaker than the forward causality, might be more sensitive to this definition than the forward causality. To assess whether the triphasic pattern is robust against the definition of the spike train order, we conducted a second analysis in which we redefined the spike train order using the Granger causality itself and evaluated the temporal pattern of the population causality. For a given cell pair, time-averaged Granger causality values (cue period) in the two directions were compared, and the larger value was defined as that in the forward direction. This procedure did not essentially change the results: the forward causality showed two separate significant peaks (supplemental Fig. 4B,E, available at www.jneurosci.org as supplemental material), and the backward causality revealed slow but significant development in the lower-frequency range (supplemental Fig. 4C,E, available at www.jneurosci.org as supplemental material). Therefore, the observed causality dynamics in both the forward and backward directions were not likely the result of an erroneous assignment of the directions.
To test whether the observed dynamics of Granger causality is specific to the cell pairs with the displaced CCG peak, we examined the dynamics of Granger causality for cell pairs that showed a center peak on the CCG (n = 14). As in the above analysis, directions of the causality for individual cell pairs were assigned on the basis of the time-averaged causality value. In contrast to the displaced peak pairs, Granger causality dynamics of the center peak pairs was quite similar between the forward and backward directions: they only showed an early component, and did not show a slowly developing component (supplemental Fig. 5B,C, available at www.jneurosci.org as supplemental material). The forward causality was slightly stronger than the backward causality, likely attributed simply to the assignment of each direction on the basis of the strength of time-averaged causality itself. The gamma Granger causality in both directions (30–65 Hz) reached significance immediately after cue onset (forward: peak, 150 ms; p < 0.003, permutation test; n = 14; backward: peak, 150 ms; p < 0.02) (supplemental Fig. 5E, available at www.jneurosci.org as supplemental material), and the directionality of the gamma Granger causality was not significant at any time point (p > 0.08) (supplemental Fig. 5F, available at www.jneurosci.org as supplemental material). Half-maximum latencies of the forward and backward causalities were not significantly different from each other (p > 0.7, paired t test). Peak frequencies of the forward and backward causality were also not significantly different from each other (p > 0.1). These results suggest that the triphasic dynamics of directional influences was specific to the cell pairs with the displaced CCG peak. Furthermore, from the above analysis, it seems unlikely that the slowly developing backward causality observed for displaced CCG peak pairs was attributable to a contribution of contaminated center peak pairs.
Recently, it has been shown that MS affects the firing rate (Leopold and Logothetis, 1998; Martinez-Conde et al., 2000; Bosman et al., 2009) and spike–LFP coherence (Bosman et al., 2009) in areas V1 and V4. This raises the possibility that the complex dynamics of Granger causality we observed might be attributable to its modulation by MSs. To assess this possibility, we detected the MSs within the fixation window during the cue period and examined whether the MSs affected the observed neuronal responses including the firing rate, coherence, and Granger causality. We conducted this analysis for all the recording sessions from one monkey whose eye-trace data were stored for the analysis (12 of 21 cell pairs with a displaced CCG peak were included).
We first calculated the MS-triggered average of the firing rate as a function of time around the MSs. To control for the effect of the stimulus-locked response, we conducted the same calculation after trial shifting of spike trains against MS events. This trial-shifted control was then subtracted to compute the shift-subtracted, MS-triggered average of the firing rate for evaluating the changes in the firing rate by MSs. We found significant changes in the firing rate that were triggered by MSs (maximum modulation: 2.8 Hz, 185 ms after MS onset; p < 0.008, permutation test) (supplemental Fig. 6A, available at www.jneurosci.org as supplemental material). This result is in line with previous observations reported in areas V1 and V4 (Leopold and Logothetis, 1998; Martinez-Conde et al., 2000; Bosman et al., 2009). To evaluate the practical impact of this effect on the actual neuronal responses during the task, we compiled the time course of MS-induced modulation of the firing rates by convolving the timing of MS onset with the shift-subtracted, MS-triggered average of the firing rate. Again, we found significant MS-induced modulation of the firing rates in the cue period (maximum modulation: −0.33 Hz, 575 ms after cue onset; p < 0.002, permutation test) (supplemental Fig. 6B, available at www.jneurosci.org as supplemental material). It should be noted, however, that these modulations were far smaller than the visual response to the presented stimulus itself (20–30 Hz) (see Fig. 5D), as expected in the IT cortex (Leopold and Logothetis, 1998).
We then examined the effects of MSs on the coherence and/or Granger causality by conducting the same analyses as those for the firing rate. Shift-subtracted, MS-triggered averages of coherence and Granger causality were found to be smaller than those signals elicited by visual stimuli (supplemental Fig. 7A–C, available at www.jneurosci.org as supplemental material). MS-induced modulation of coherence and Granger causality in the cue period (supplemental Fig. 7D–F, available at www.jneurosci.org as supplemental material) did not reveal prominent components at any time point after cue onset, in contrast to those that were elicited by visual stimuli (compare supplemental Fig. 7D–F, available at www.jneurosci.org as supplemental material, with Fig. 4F–H). We next calculated the gamma components of the shift-subtracted, MS-triggered averages and the MS-induced modulations for both coherence and Granger causality (supplemental Fig. 8, available at www.jneurosci.org as supplemental material). The gamma components of the shift-subtracted, MS-triggered average of coherence and Granger causality were not significant at any time point from MS onset (coherence, p > 0.2; forward Granger causality, p > 0.9; backward Granger causality, p > 0.06, permutation test) (supplemental Fig. 8A,B, available at www.jneurosci.org as supplemental material), and their maximum amplitudes were significantly smaller than those that were elicited by visual stimuli (coherence, p < 0.001; forward Granger causality, p < 0.004; backward Granger causality, p < 0.001, two sample t test) (compare supplemental Fig. 8A,B, available at www.jneurosci.org as supplemental material, with Fig. 5A,B). In the cue period, the gamma components of the MS-induced modulation of coherence and Granger causality were not significant at any time point after cue onset (coherence, p > 0.1; forward Granger causality, p > 0.5; backward Granger causality, p > 0.06, permutation test) (supplemental Fig. 8C,D, available at www.jneurosci.org as supplemental material), and their maximum amplitudes were, again, significantly smaller than those that were elicited by visual stimuli (p < 0.001 for coherence and forward/backward Granger causality, two-sample t test) (compare supplemental Fig. 8C,D, available at www.jneurosci.org as supplemental material, with Fig. 5A,B). Note that neither the MS-induced modulation of coherence or Granger causality in the gamma range revealed even a small temporal pattern such as we observed for those that were elicited by visual stimuli. Therefore, we conclude that MSs cannot explain the stimulus-elicited temporal pattern of gamma coherence or Granger causality between SUAs in the IT cortex.
Stimulus-dependence of the Granger causality dynamics
The above data showed a temporally dynamic feature of the Granger causality between spikes from single neurons. The next question is whether these dynamics of the directional influences occur in a stimulus-dependent manner. In a previous study, we prepared 64,000 FOs and their corresponding NFOs and demonstrated that spike-timing correlation between pairs of neurons was significantly stronger for their pFOs compared with the corresponding pNFOs composed of the same set of parts. The firing rates were not significantly different between pFOs and pNFOs (Hirabayashi and Miyashita, 2005). We thus examined whether the dynamics of Granger causality also depends on the feature configuration within a whole object. As a population, pFOs elicited strong forward causality in the early and late epochs of the cue period and slowly developing backward causality (Fig. 7, left). In contrast, Granger causality elicited by the corresponding pNFOs was much weaker for the same cell pairs and only transiently appeared in the early epoch, in the forward direction (Fig. 7, right).
Stimulus dependence of Granger causality dynamics. Each panel represents population average (n = 21) of Granger causality dynamics in each direction in response to optimal FOs or NFOs. Optimal FO and NFO for a given cell pair were composed of the same set of facial feature-like parts that were arranged in a face-like pattern and in a random pattern, respectively.
Quantitative analysis of gamma Granger causality dynamics showed that pFOs elicited temporally separated two significant components of forward causality (30–80 Hz; peak: 150 ms and 525 ms; p < 0.002 and p < 0.02, respectively, permutation test) (Fig. 8A, left, red trace) and elicited a single, slowly developing significant component of backward causality (30–65 Hz; peak: 475 ms; p < 0.005) (Fig. 8A, left, blue trace). On the contrary, pNFOs elicited only an early transient forward causality (peak: 75 ms; p < 0.03) (Fig. 8A, right, red trace). No systematic difference was observed between the population firing rates for pFOs and pNFOs at any time point (Fig. 8B), suggesting that the observed dependence of Granger causality dynamics on the presented stimulus was not attributable to the firing rate difference. We conducted two-way ANOVA to quantitatively investigate the difference in the Granger causality between pFOs and pNFOs (Fig. 8C). The pFO-evoked forward causality was maximum at 150 ms (from cue onset) in the early epoch when the pFO-evoked forward causality reached significance. At that time point, pFO-evoked forward causality was significantly greater than that for pNFOs or than the backward component of the pFO-evoked causality (p < 0.04, two-way ANOVA followed by LSD) (Fig. 8C, left). The pFO-evoked backward causality was the maximum at 475 ms in the middle epoch when the backward component of the pFO-evoked causality reached significance. At that time point, the forward component of pFO-evoked causality was significantly larger than that for pNFOs (p < 0.01), whereas the backward component was not significantly different between pFOs and pNFOs (p > 0.2) (Fig. 8C, middle). The pFO-evoked forward causality was maximum at 525 ms in the late epoch when the pFO-evoked forward causality again reached significance. At that time point, pFO-derived forward causality was significantly larger than that for pNFOs or than the backward component of the pFO-evoked causality (p < 0.01) (Fig. 8C, right), as was observed in the early epoch. In contrast to the stimulus dependence of the Granger causality, there was no significant difference between the firing rates in response to pFOs and pNFOs at the above time points (pFO vs pNFO: early, 18.3 ± 1.4 vs 16.5 ± 1.9 Hz, mean ± SEM; p > 0.06, paired t test; middle, 19.5 ± 1.7 vs 18.4 ± 1.8 Hz; p > 0.4; late, 18.8 ± 1.7 vs 17.7 ± 1.8 Hz; p > 0.4). Together, these results indicate that the triphasic dynamics of directional coupling occurs only for pFOs. Thus, it is tempting to suggest that stimulus-dependent unidirectional causality in the early and late epochs might reflect the signal processing regarding different aspects of the object recognition, such as structural and semantic analyses.
Statistical analysis of stimulus-dependent Granger causality dynamics. A, Population average time courses of gamma Granger causality elicited by optimal FOs and NFOs (n = 21). Thin traces are mean ± SEM. Red and blue horizontal lines, Confidence limit for forward and backward directions, respectively (p = 0.05, corrected for multiple comparisons). B, Dynamics of population mean firing rate for optimal FOs and NFOs. Thin traces are mean ± SEM. C, Population averages of gamma Granger causality at the time point when FO-evoked causality was the maximum in each epoch. Each epoch was defined as that in which FO-derived gamma Granger causality in either direction was significant (see Results). Error bars represent SEM. Asterisks indicate statistical significance assessed by two-way ANOVA followed by LSD.
Discussion
In the present study, we examined the Granger causality between spike trains of single neurons in macaque IT cortex that exhibited a significant displaced peak on the CCG, suggesting pauci-synaptic serial connections in a classical view. In response to the stimulus, these pairs of neurons showed Granger causality in the gamma range, whose directionality was consistent with that of the CCG peak. Temporal dynamics of the Granger causality exhibited a triphasic pattern: a transient forward component was followed by a slowly developing backward component and subsequent reappearance of the forward component. These temporal dynamics of Granger causality depended on the feature configuration within a presented object and were not observed for cell pairs that exhibited a center peak on the CCG. This multiphasic flow of information is likely to be a reflection of the multistage network processing in the recognition of visual objects in the IT cortex (Miyashita, 2004).
Granger causality analysis for spike trains of single neurons
To apply the parametric Granger causality analysis (Ding et al., 2006) to point processes such as spike trains, they have to be converted into continuous data for conducting MVAR modeling of the signal. This conversion might lead to loss of information as a point process inherent in the original spike trains (Truccolo et al., 2005; Dhamala et al., 2008a,b; Nedungadi et al., 2009). Nonparametric Granger causality analysis bypasses the MVAR modeling of the signal (Dhamala et al., 2008a,b) and thus can be applied directly to spike trains (Nedungadi et al., 2009). In previous studies, it has been confirmed that nonparametric Granger causality analysis for LFP and simulated data yielded qualitatively similar results as those obtained by parametric Granger causality analysis (Dhamala et al., 2008a,b). With respect to spike trains, it has been confirmed that nonparametric Granger causality analysis correctly estimated the anatomically identified unidirectional connectivity (Nedungadi et al., 2009). In the present study, for both the parametric and nonparametric methods, cell pairs with a displaced CCG peak showed a significantly forward-dominant Granger causality in the cue period, and trial shifting significantly eliminated this causality. Thus, with the aforementioned previous studies, we suggest that our results do not depend on the used method for estimating the Granger causality.
In the analysis of neuronal interactions including cross-correlation, coherence, and Granger causality, stationarity of the signal is a prerequisite for their correct estimation (Ding et al., 2000; Jarvis and Mitra, 2001; Brown et al., 2004; Young and Eggermont, 2009). Analyzing the interactions in a sliding short time window, in which the signal could be treated as locally stationary, is an effective and practical strategy to cope with this issue (Ding et al., 2000; Jarvis and Mitra, 2001; Brown et al., 2004). In the present study, we also calculated the Granger causality using a sliding short time window, with the length of ±125 ms, which falls within the range used in previous studies of spectral analysis for spikes or LFPs (Fries et al., 2008; Zhou et al., 2008). We also calculated the Granger causality using trial-shifted data to determine whether the calculated Granger causality simply resulted from the firing dynamics of each neuron. Trial shifting almost completely diminished the observed Granger causality and its directionality (Figs. 1, 2). All the calculations of the Granger causality analysis in the present study were evaluated as a comparison between the results using raw and trial-shifted data, or were evaluated after subtraction of the trial-shifted results. These procedures eliminated the possibility that the observed Granger causality was attributable to the nonstationarity of the signal, including a sequential increase and/or decrease in the firing rates without connectivity.
What can be observable as the effect of the Granger causality on the spike trains? Kayser and Logothetis (2009) demonstrated that the strength of LFP–LFP Granger causality from the superior temporal sulcus to the auditory cortex was positively correlated with the amplitude of multiunit activity (MUA) modulation in the auditory cortex during a combined auditory and visual stimulation. Pesaran et al. (2008) examined the spike–LFP coherence between the dorsal premortor area (PMd) and the parietal reach region (PRR) during a cue-triggered decision-making task. The results demonstrated that the firing rates at a pair of recording sites that revealed significant coherence indicated the animal's decision earlier than the firing rates at another pair of recording sites that did not exhibit significant coherence. These results suggest that one of the observable effects of the Granger causality is firing rate modulation of the neuron that is the target of the causality. In the present study, however, we did not observe such modulation of the firing rate that was time locked to the transient Granger causality (data not shown). This may have occurred, at least in part, because the firing rate of a single neuron reflects inputs from a large number of other neurons, which could cause difficulties in extracting only the effect of the inputs from a specific single neuron. This might be different from a case in which the Granger causality and its effect are examined between MUAs or LFPs. It has been reported that the strength of coherence tends to be weaker when the analyzed signals are SUAs compared with LFPs or MUAs (Zeitler et al., 2006; Gregoriou et al., 2009). Considering that individual neurons act as components in a network in which they are embedded, even a weak impact of each interaction between single neurons can affect the network behavior as a population.
Coherence and Granger causality
Coherence is a widely used measure to assess the interrelationships between neuronal signals. Individual spike–spike coherences in the present study revealed distinct peaks with various peak frequencies in the gamma range. Several studies have reported variance of the peak frequency for individual gamma coherence. Jutras et al. (2009) reported that the frequency of hippocampal SUA–LFP coherence is separated into low and high gamma ranges. Sirota et al. (2008) also reported that the peak frequencies of neocortical SUA–SUA and SUA–LFP coherence are distributed over a wide range within the gamma band.
The frequency of gamma oscillation depends on the surrounding circuitry and its behavior (Atallah and Scanziani, 2009; Mann and Mody, 2010). Therefore, pair-by-pair variability of the peak frequencies might be expected to be larger in SUA–SUA coherence than those between MUAs or between MUA and LFP, since the variability of the surrounding circuit would be larger in the case of single-unit pairs compared with other configurations.
The phase of coherence represents the time lag between the signals and thus has been used as an index of directionality of the interaction (Kocsis and Vertes, 1997; Bressler and Kelso, 2001; Saalmann et al., 2007). In the present study, the directionality index of Granger causality was significantly correlated with the AI of the CCG or with the phase of coherence (Fig. 3C–E and supplemental Figs. 2, 3, available at www.jneurosci.org as supplemental material). These data suggest that the obtained results of nonparametric Granger causality analysis were consistent with the results of both cross-correlation and coherence analyses in terms of the directionality, further supporting the reliability of the present analysis of Granger causality.
In the present study, the dynamics of Granger causality was more rapid and complex than that of the coherence: whereas the gamma coherence remained to be significant throughout the cue period, the gamma Granger causality and its directionality revealed more dynamic changes (Figs. 4, 5). The amplitude of coherence provides the strength of interrelationship regardless of its directionality. Thus, the dynamics of neuronal interactions obtained by coherence analysis might, in some cases, be further decomposed into more rapid and complex, directed interactions as shown in the present study. It should be noted, however, that spike–LFP coherence can be used as another method to examine the directed interactions between distant locations, considering that spike and LFP correspond to the signals related to neuronal output and input, respectively. Pesaran et al. (2008) examined the interactions between the areas PMd and PRR using this technique. They found that the PMd spike–PRR LFP coherence appeared first after cue onset and the opposite PMd LFP–PRR spike coherence appeared later, suggesting that information transferred first from PMd to PRR and then back from PRR to PMd.
Cross-correlation and Granger causality
In a classical view of functional connectivity, a displaced CCG peak suggests pauci-synaptic serial connections, whereas a center peak on the CCG suggests common inputs or recurrent connections (Perkel et al., 1967; Moore et al., 1970). However, cortical neurons are embedded in a complex network, and thus the impacts of surrounding local circuit also have to be considered to explain the interactions between the recorded pair of neurons (Zhou et al., 2008; Ostojic et al., 2009). Indeed, the present triphasic causality dynamics observed for cell pairs with a displaced CCG peak suggests a nonclassical view of the functional connectivity within a cortical local circuit.
From a technical point of view, however, it is difficult to distinguish even the two classical types of functional connectivities using a CCG. For example, a given pair of neurons in a local cortical circuit often involves both types of interactions (Yoshimura et al., 2005). Additionally, even common inputs or recurrent connections might yield a displaced CCG peak if the conduction time or the connection strength was different between the recorded pair of cells. These possibilities raise the concern that the observed time-frequency pattern of Granger causality estimated for displaced peak pairs might result from non-negligible direct contributions of common inputs or recurrent connections. However, when we analyzed only the neuron pairs whose CCG peaks were classified as the center peak, the time-frequency pattern of Granger causality was largely different from that observed for cell pairs with a displaced CCG peak (supplemental Fig. 5, available at www.jneurosci.org as supplemental material) (i.e., forward and backward components with similar frequency range appeared after cue onset with a short latency, and both the components simply decayed monotonically without a noticeable rebound). These patterns do not explain the late reappearance of forward causality nor the delayed development of backward causality observed for displaced CCG peak pairs. In addition, for the displaced CCG peak pairs, the time lag between the early component of the forward causality and the slowly developing backward causality was over 100 ms (Figs. 4, 5), which should be far beyond the difference in the conduction delay of common inputs to the recorded pair of cells. Furthermore, we also confirmed that the slowly developing backward causality remained significant even if we only analyzed the cell pairs whose CCG showed the AI of 1.0, the upper limit of the index (i.e., most displaced). We thus conclude that it is unlikely that the observed dynamics of directional interdependence for displaced CCG peak pairs were merely the result of direct contributions of simple common inputs or recurrent connections.
Is the backward causality detectable in the CCG? In supplemental Fig. 9 (available at www.jneurosci.org as supplemental material), we demonstrate an example of the backward causality reflected on the CCG peak for a cell pair shown in Figure 1. This cell pair exhibited unidirectional forward causality at 225 ms after cue onset and showed a slowly developing backward causality at 475 ms, as typically observed (supplemental Fig. 9, insets, available at www.jneurosci.org as supplemental material). For this cell pair, we compiled corresponding CCGs using the spike trains in the same time windows. At 225 ms after cue onset, a primary displaced peak was observed, with the directionality expected from the Granger causality (supplemental Fig. 9A, available at www.jneurosci.org as supplemental material). At 475 ms, an additional secondary peak appeared on the opposite side of the CCG, which appears to correspond to the development of backward causality (supplemental Fig. 9B, available at www.jneurosci.org as supplemental material). These results indicate that backward causality can be also observed on the CCG, in some cases.
In most cases, however, a sliding-window CCG leads to a shortage of spikes because of the shortening of the time window, which reduces the reliability and detectability of the CCG analysis. Furthermore, the detectability of neuronal interactions is, in general, improved in the spectral analysis compared with the conventional cross-correlation analysis (Jarvis and Mitra, 2001; Fries et al., 2008). Therefore, the present Granger causality analysis will provide a new and powerful approach to uncover the rapid dynamics of directed interactions beyond a classical view of functional connectivity inferred from a CCG. Especially, the dynamics of directed information flow between single neurons whose spiking activities are characterized in a cognitive task would be an intriguing issue for future study.
Possible mechanisms underlying the directed interactions
Differential dynamics of the forward and backward causality in the present results (Figs. 4, 5) implies that different neuronal circuits are involved in these interactions. In the present study, we analyzed the directed interactions between neurons that exhibited a significant displaced peak on the CCG. A straightforward interpretation regarding the mechanisms underlying the observed forward causality is therefore pauci-synaptic unidirectional interactions between the neurons. It has been shown previously that nonparametric Granger causality analysis faithfully represented anatomically identified thalamo-cortical monosynaptic connections (Nedungadi et al., 2009). This interpretation is also supported by the observation that the forward causality appeared in the early phase of the response, even earlier than the timing of the peak firing rate (Figs. 4, 5).
The mechanism underlying the observed backward causality, however, is likely to be substantially different. The latency of the backward causality was significantly longer than that of the forward causality, suggesting that the path length underlying this backward causality is significantly longer than that underlying the forward causality. Zhou et al. (2008) showed that stimulus-evoked gamma coherence between SUAs in V1 developed slowly over 500 ms after cue onset, and discussed that this slow development of the coherence is likely attributable to the network activity of surrounding neurons. Samonds et al. (2009) demonstrated that spike-timing correlation in V1 showed temporally discrete early and late increase in response to the visual stimulus, and discussed that the circuits underlying these early and late increase in the neuronal interactions might be local and distant, respectively. Therefore, the observed backward causality might be the result of the interactions between the recorded neurons and their surrounding network of neurons responding to the same stimulus. The delayed backward causality suggests that the local cortical circuit includes a closed-loop structure, an intermediate between simple feedforward and recurrent configurations, that might enable various computations (Ahissar and Kleinfeld, 2003). We do not, however, exclude the possibility that more distant circuitry is involved in the observed backward causality, such as top-down signaling from other cortical areas (Tomita et al., 1999; Naya et al., 2001; Miyashita, 2004; Takeda et al., 2005). Feedback activity, which can be reflected in the backward causality, might play a critical role in enhancing the transmission of information as shown in the thalamocortical pathway (Briggs and Usrey, 2008).
Functional implications of the temporal pattern of information flow
Temporally discrete patterns of activity and their functional implications have been reported for other neuronal signals. In human event-related potential (ERP) studies, it has long been considered that ERP signals with different latencies in the visual cortex code different information about the presented visual stimulus. In visual word recognition, for example, temporally distinct components have been associated with different processing stages including orthographic and lexical processing (Holcomb and Grainger, 2006; Patterson et al., 2006). In face recognition, it has been hypothesized that structural information is processed first, followed by the processing of semantic information (e.g., familiarity or identification) (Bruce and Young, 1986; Tovee and Cohen-Tovee, 1993). Indeed, several lines of evidence have been shown regarding differences in the information content between early and late components of the ERP signal: an early component at ∼170 ms is thought to reflect structural processing, whereas a late component at ∼400–600 ms reflects semantic processing (Bentin and Deouell, 2000; Eimer, 2000). In the present study, we found three major components of Granger causality: an early forward (peak, 150 ms) component, a middle backward (400 ms) component, and a late forward (550 ms) component (Fig. 5). The early and late forward components were predominantly observed in response to FOs (Fig. 8), implying that they might be associated with the structural and semantic processing (and/or with the perceptual decision-making regarding these information) of the presented FO, respectively. Note, however, the observed triphasic dynamics might be one of possible patterns of more general multiphasic dynamics of directed interactions, which might depend on the behavioral conditions including the cognitive task or the presented stimulus.
The strength of coupling between single neurons in the visual cortex has been shown to be modulated depending on the presented stimulus or behavioral context (Engel et al., 2001; Salinas and Sejnowski, 2001; Hirabayashi and Miyashita, 2005; Kohn and Smith, 2005). It has not been examined, however, whether the directionality of information flow between single neuron changes in a stimulus-dependent manner. The present results provide the first example of the Granger causality dynamics between the spike trains of single neurons in behaving monkeys. Further elucidation of the mechanisms and functions of this rapid and complex dynamics of directional interdependence will be an intriguing issue for future study.
Footnotes
This work was supported by a Grant-in-Aid for Specially Promoted Research (19002010) to Y.M. from the Ministry of Education, Culture, Sports, Science, and Technology of Japan and by the Takeda Science Foundation. We thank Dr. G. Rangarajan for kindly providing the script to conduct the nonparametric Granger causality analysis, Drs. Y. Sakurai and S. Takahashi for kindly providing the script of ICSort, and Dr. M. Takeda for helpful discussion.
- Correspondence should be addressed to Dr. Yasushi Miyashita, Department of Physiology, The University of Tokyo School of Medicine, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-0033, Japan. yasushi_miyashita{at}m.u-tokyo.ac.jp
