A convolutional neural network for estimating synaptic connectivity from spike trains

The recent increase in reliable, simultaneous high channel count extracellular recordings is exciting for physiologists and theoreticians because it offers the possibility of reconstructing the underlying neuronal circuits. We recently presented a method of inferring this circuit connectivity from neuronal spike trains by applying the generalized linear model to cross-correlograms. Although the algorithm can do a good job of circuit reconstruction, the parameters need to be carefully tuned for each individual dataset. Here we present another method using a Convolutional Neural Network for Estimating synaptic Connectivity from spike trains. After adaptation to huge amounts of simulated data, this method robustly captures the specific feature of monosynaptic impact in a noisy cross-correlogram. There are no user-adjustable parameters. With this new method, we have constructed diagrams of neuronal circuits recorded in several cortical areas of monkeys.


Results
Training and validating with synthetic data. CoNNECT infers the presence or absence of monosynaptic connections between a pair of neurons and estimates the amplitude of the postsynaptic potential (PSP) that one neuron would drive in another. The estimation is performed by applying a convolutional neural network [15][16][17][18] to a cross-correlogram obtained for every pair of spike trains (Fig. 1a). The network has an output layer consisting of two units. One unit indicates the presence or absence of connectivity with a real value z ∈ [0, 1] by thresholding at 0.5. Another is the level of PSP represented in a unit of (mV). The network was trained with spike trains generated by a numerical simulation of a network of multiple-timescale adaptive threshold (MAT) model neurons [21][22][23][24][25] interacting through fixed synapses. In a large-scale simulation, we applied fluctuating inputs to a subset of neurons to reproduce large fluctuations in real spike trains in vivo (Fig. 1b). Figure 1c,d, and e demonstrate sample spike trains, histograms of the firing rates of excitatory and inhibitory neurons 26 , and firing irregularity measured in terms of the local variation of the interspike intervals Lv 27,28 . The training data does not contain many low firing rate neurons, considering the situation that low firing units are often discarded when analyzing real data. The details of the learning procedure are summarized in "Methods" section.
We validated the estimation performance of CoNNECT using novel spike trains generated by another neuronal circuit with different connections. Figure 2a depicts an estimated connection matrix, referenced to the true connection matrix, of 50 neurons. Here, the estimation was done with spike trains recorded for 120 min. Of 50 spike trains, 40 and 10 are, respectively, sampled from 800 excitatory and 200 inhibitory neurons. Figure 2b compares the estimated PSPs against true values. We have presented an estimated PSP as being 0 if the connection is not detected. Points lying on the nonzero x-axis are existing connections that were not detected or FNs. Points lying on the nonzero y-axis are spurious connections assigned for unconnected pairs or FPs. Figure 2c depicts how the numbers of FNs and FPs for excitatory and inhibitory categories changed with the recording duration or the length of spike trains (10,30, and 120 min). While the number of FPs or spurious connections does not depend largely on the recording duration, the number of FNs or missing connections decreased with the period, implying that more synaptic connections of weaker strength are revealed by increasing the recording time.
Comparison with other estimation methods. There are many algorithms that were developed to estimate synaptic connections from spike trains [1][2][3][4][5][6][7][8][9][10][11][12][13]19,[29][30][31] . We compared CoNNECT with the conventional cross-correlation method (CC) 19 , Jittering method 4 , Extended GLM 13 , and GLMCC 21 for their ability to estimate connectivity using synthetic data. Figure 3 shows connection matrices determined by the four methods referenced to the true connection matrices. In the lower panels, we demonstrated the performances in terms of the false positive (discovery) rate (FPR) and false negative (omission) rate (FNR) for excitatory and inhibitory categories; smaller values are better. Here we estimated the mean and SD of the performance by applying each method to 8 test datasets of 50 neurons. Overall performance with FPR and FNR was measured in terms of the Matthews correlation coefficient (MCC) (see "Methods" section). The MCCs for these estimation methods are shown in the right edge panel; a larger MCC is better. For evaluating the performances, we adopted spiking data generated by a network of MAT models and a network of Hodgkin-Huxley (HH) type models ("Methods" section). In computing the numbers of FPs and FNs, we ignored small excitatory connections, which are inherently difficult to discern with this observation duration. We took the lower thresholds as 0.1 mV for the MAT simulation and 1 mV for the HH simulation so that the visible connectivity is about 10%, but the relative performances between different models are unchanged even if we change the thresholds.
The conventional cross-correlation analysis produced many FPs, revealing a vulnerability to fluctuations in cross-correlograms. The Jittering method succeeded in avoiding FPs but missed many existing connections, thus generating many FNs. The Extended GLM method of given parameters was also rather conservative. In comparison to these methods, GLMCC and CoNNECT have better performance, producing a small number of FPs and FNs and a larger MCC value. Here we have modified GLMCC so that it achieves higher performance than the original algorithm 14 by using the likelihood ratio test to determine the statistical significance ("Methods" section). When comparing these two algorithms, GLMCC was slightly conservative, producing more FNs, while CoNNECT tended to suggest more connections, producing more FPs.
The converted GLMCC was better than CoNNECT for the HH model data (Fig. 3b), but the converse was true for the MAT model data (Fig. 3a). This might be because CoNNECT was trained using the MAT model data of a similar kind, and GLMCC was constructed by considering the HH model simulation. Although the model Comparison of different learning conditions. While the convolutional neural network has the advantage that tens of thousands of parameters can be suitably adjusted to reproduce given datasets, it does not guarantee the generalization capability. We evaluated the generalization capability of our convolutional network by changing the number of out-channels representing the degree of system adaptability from 1 to 10. Figure 4a depicts the numbers of FPs and FNs in the above, and the overall performances measured in terms of MCC in the below, which were obtained for the MAT model simulation data (left panel) and the Hodgkin-Huxley model simulation data (right panel). We observe that the convolutional network consisting of 1 channel slightly "under-fits" because of the little flexibility, whereas that of 10 channels slightly over-fits the data, exhibiting slightly lower MCC. Thus we have employed the network of 5 channels, consisting of about fifty thousand parameters.   To make the convolutional network applicable to data of a wider variety, we have trained the network using the cross-correlograms augmented by rescaling the time. Figure 4b depicts the estimation performances of networks trained using the cross-correlations rescaled by 1/4 and 1/2 times of the original (indicated as "1/4" and "1/2"), the original cross-correlations ("1"), and all the data ("1/4 + 1/2 + 1"). The networks trained with lower firing rates exhibited lower performances. We have adopted the network trained with all the data ("1/4 + 1/2 + 1") because it gave the highest performance in estimating connectivity.

Cross-correlograms.
To observe the situations in which different estimation methods succeeded or failed in detecting the presence or absence of synaptic connectivity, we examined sample cross-correlograms of neuron pairs of a network of MAT model neurons.  www.nature.com/scientificreports/ in real biological data. These were produced by external fluctuations added to a subset of neurons, making the connectivity inference difficult. The inference results obtained by the four estimation methods are distinguished with colors; magenta, cyan, and gray represent that estimated connections were excitatory, inhibitory, or unconnected, respectively. We also superimposed a GLM function fitted to each cross-correlogram. Figure 5a,b depict sample cross-correlograms of neuron pairs that are connected with excitatory and inhibitory synapses, respectively. For the first three cross-correlograms from the top, all four estimation methods succeeded in detecting excitatory or inhibitory connections, thus making true positive (TP) estimations. For the fourth case, the Jittering method failed to detect the connection. This implies that the Jittering method is rather conservative for producing FPs, and as a result, has produced many FNs. In this case, the cross-correlation method (CC) has mistaken the excitatory synapse as inhibition due to the large wavy fluctuation in the crosscorrelogram. For the last cases, all four estimation methods failed to detect the connection, resulting in FNs. This would have been because the original connections were not strong enough to produce significant impacts on the cross-correlograms.     Preprocessing experimental data. Some of the experimentally available cross-correlograms exhibit a sharp drop near the origin for a few ms due to the shadowing effect, in which near-synchronous spikes cannot be detected 32 . This effect disrupts the estimation of synaptic impacts that should appear near the origin of the crosscorrelogram. The data were obtained with a sorting algorithm specifically used for the Utah array exhibit rather broad shadowing effects larger than 1 ms (up to 1.75 ms). Here, we analyzed the experimental data by removing an interval of 0 ± 2 ms in the cross-correlogram and applying the estimation method to a cross-correlogram obtained by concatenating the remaining left and right parts (Fig. 6a,b). We also conducted this operation in the analysis of synthetic data. Figure 6c demonstrates the cross-correlograms of sample neuron pairs for which both CoNNECT and GLMCC estimated connections which were excitatory, inhibitory, or absent (unconnected). It was observed that the real cross-correlograms are accompanied by large fluctuations. Nevertheless, CoNNECT and GLMCC are able to detect the likely presence or absence of synaptic interaction by ignoring the severe fluctuations.   , were applied to cross-correlograms. Their estimation (excitatory, inhibitory, and unconnected) are respectively distinguished with colors (magenta, cyan, and gray). The lines plotted on the cross-correlograms are the GLM functions fitted by GLMCC. The causal impact from a pre-neuron to a post-neuron appears on the right half in each cross-correlogram.  Figure 7 depicts the estimated connections for the entire three datasets of PF, IT, and V1. The units in the connection matrices are arranged in the order provided by a sorting algorithm, and accordingly, units of neighboring indexes of the matrices tended to have been spatially closely located. All three connection matrices had more components in near diagonal elements, implying that neurons in a nearby location are more likely to be connected. The firing rate and irregularity (the local variation of the interspike intervals Lv 27,28 ) are shown in the rightmost panels. The summary statistics in Table 1 reflect differences in firing rate between excitatory and inhibitory cells in PF and IT but not V1. The firing irregularity of excitatory neurons is slightly higher than that of inhibitory neurons, consistent with the previous results. Table 1 summarizes the statistics of the three datasets. Each neuron is assigned as putative excitatory, putative inhibitory, or undetermined, according to whether the excitatory-inhibitory (E-I) dominance index is positive, negative, or undetermined (or zero), respectively. Here, the E-I dominance index is defined as  www.nature.com/scientificreports/ few, the majority of d ei is either 0 or 1. Though we have obtained many connections, the total number of all pairs is enormous, scaling with the square of the number of units, and accordingly, the connectivity is sparse (less than 1% for each (directed) pair of neurons). In contrast to synthetic data, the currently available experimental data do not contain information regarding the true connectivity. To examine the stability of the estimation, we split the recordings in half and compared estimated connections from each half. If the real connectivity is stable, we may expect the estimated connections have overlap between the first and second halves. Figure 8a represents the connection matrices obtained from the first and second halves of the spike trains recorded from PF, IT, and V1. Figure 8b compares the estimated PSPs in two periods. Many estimated connections appear only on one of the two. This might be simply due to statistical fluctuation or due to real changes in synaptic connectivity. Nevertheless, it may be noteworthy that the excitatory connections of large amplitudes were detected relatively consistently between the first and second halves. Namely, they appear in the first and third quadrants diagonally, implying that they have the same signs with similar amplitudes.

Discussion
Here we have devised a new method for estimating synaptic connections based on a convolutional neural network. While this method does not require adjusting parameters for individual data, it robustly provides a reasonable estimate of synaptic connections by tolerating large fluctuations in the data. This high performance was obtained by training a convolutional neural network using a considerable amount of training data generated by simulating a network of model spiking neurons subject to fluctuating current. We compared CoNNECT with the conventional cross-correlation method, the Jittering method, Extended GLM, and GLMCC in their ability to estimate connectivity, using synthetic data obtained by simulating neuronal circuitries of fixed synaptic connections. Both CoNNECT and GLMCC exhibited high performance in predicting individual synaptic connections, superior to other methods.
Then we applied CoNNECT to simultaneously recorded spike trains recorded from monkeys using the Utah arrays. We have found that the connections among recorded units are sparse; they are less than 1% for all three datasets. To test the reliability of the estimation, we divided the entire recording interval in half and estimated connections for respective intervals. We have seen that strong excitatory connections overlap between the periods. This result implies that the estimation is reliable for the strong connectivity, and the connectivity lasts at least for hours.
The cross-correlograms of real biological data (Fig. 6) turned out to be even more complicated than those of synthetic data (Fig. 5), which were generated by adding large fluctuations to individual neurons (Fig. 1). The complicated features in real cross-correlograms were not only due to fluctuations in real circuitry but also due to the sorting algorithm. The most severe bottleneck in estimating connectivity may have been the shadowing effect of a few ms, in which near-synchronous spikes were not detected (Fig. 6a); this effect might hide the first part of a monosynaptic impact, which is expected to show up in a few ms in a cross-correlogram. If the sorting algorithm is improved such that the shadowing duration is shortened, the estimation might be more reliable.
In this study, we have employed the convolutional neural network to capture the specific signature of monosynaptic impact in a cross-correlogram image. While the convolutional network is known to be robust against the translation of images, the monosynaptic impact is expected to appear at a specific location in the cross-correlogram, particularly exhibiting the delay of a few milliseconds. Thus it might be an interesting challenge to search for other learning algorithms that utilize such information and perform better than the convolutional network.
We used data augmentation technique 33 to increase the number of training examples artificially. Data augmentation is known to improve the performance on various tasks in computer vision 18,34 and acoustic signal processing 35,36 . Here we augmented the cross-correlation data by rescaling the time to capture the diverse synaptic interactions. This augmentation also improved the estimation performance of synaptic connectivity (Fig. 4b).
Recently, several authors proposed a more systematic approach for data augmentation, e.g., generating augmented data using generative adversarial networks (GANs) 37,38 and learning the data augmentation policy 39 . Though these approaches focus on the image classification task and require a vast computational resource, it may be interesting to apply these techniques to pursue an advanced data augmentation method for synaptic connectivity estimation.
So far, we have little knowledge about neuronal circuitry in the brain. By collecting more data from high channel count recordings and applying these reliable analysis methods to them, we shall be able to obtain information about neuronal circuitry in different brain regions and learn about network characteristics and the information flow in each area. Ultimately, we expect that we will characterize the network characteristics of different brain regions processing various kinds of information.

Configuration of a neural network for estimating synaptic connectivity.
Here we describe the details of a four-layered convolutional neural network [15][16][17][18] applied to cross-correlograms obtained for every pair of spike trains to estimate the presence or absence of a connection and its postsynaptic potential (PSP) (Fig. 1). The neural network learns to find a bump or dent in the cross-correlogram caused by a monosynaptic connection.
In particular, the input consists of 100 integer values of the spike counts in a cross-correlation histogram in an interval of [−50, 50] ms with 1 ms bin size. The network comprises a 1-dimensional convolution layer, the average pooling, and the internal layer of fully connected 100 nodes. The output layer consists of two units; one indicates the presence or absence of connectivity with a real value z ∈ [0, 1] . Another is the level of PSP represented in a unit of (mV).
Training the convolutional neural network. We ran a numerical simulation of a network of 1000 neurons interacting through fixed synapses in various conditions and trained the neural network with spike trains from 400 units selected from the entire network. Thus, we constructed cross-correlograms of about 80, 000 pairs, each assigned with the teaching signals consisting of the true information about the presence or absence of connectivity (respectively represented as z = 1 or 0) and its PSP value in either direction. The training was performed using an algorithm for first-order gradient-based optimization of stochastic objective functions, based on adaptive estimates of lower-order moments, named Adam 40 . The parameters adopted in the learning are summarized in Table 2. Details of the architecture are summarized in Table 3. Figure 9 demonstrates a set of convolutional kernels that were learned with training data. From the set of learned kernels, we can see some specific features of monosynaptic impact of a few milliseconds appearing in a cross-correlogram. It is also interesting to see a kernel exhibiting a roughly monotonic gradient (the second panel from the left). This might have worked for detrending the large slow fluctuations in the cross correlogram produced by our simulation, which aimed at reproducing real situations. www.nature.com/scientificreports/ Data augmentation. To make the estimation method applicable to data of a wider range, we performed data augmentation 18,34-37 (see 33 for review). Namely, we augmented the data by rescaling the cross-correlations by 2 and 4 times and used all the data, including the original data in the learning.

Web-application program.
A ready-to-use version of the web application, the source code, and example data sets are available at our website, https://s-shino moto. com/ CONNE CT/ and are also hosted publicly on Github, accessible via https:// github. com/ shige rushi nomoto. The simulation code is also available at this Github.

Improvement of GLMCC. Original framework of GLMCC.
In the previous study 14 Table 3. Architecture of the convolutional neural network.

Convolution layer
Kernel size 10 Number of out-channels 5 Slide range 1

Number of nodes 100
Activation function ReLU Likelihood ratio test. The likely presence of the connectivity can be determined by disproving the null hypothesis that a connection is absent. In the original model, this was performed by thresholding the estimated parameters with |Ĵ ij | > 1.57z α (τ T pre post ) −1/2 , where z α , T, pre , and post are a threshold for the normal distribution, recording time, firing rates of pre-and post-synaptic neurons. But we realized that this thresholding method might induce a large asymmetry in detectability between excitatory and inhibitory connections. Instead of a simple thresholding, here we introduce the likelihood ratio test that is a general method for testing hypothesis (Chapter 11 of 41 , see also 42 ): We compute the likelihood ratio between the presence of the connectivity J ij =Ĵ ij and the absence of connectivity J ij = 0 or its logarithm: where L * J ij = c in each case is the likelihood obtained by optimizing all the other parameters with the constraint of J ij = c . It was proven that 2D obey the χ 2 distribution in a large sample limit (Wilks' theorem) 43 . Accordingly, we may reject the null hypothesis if 2D > z α , where z α is the threshold of χ 2 distribution of a significance level α . Here we have adopted α = 10 −4 .
Model validation. The performance of CoNNECT was evaluated using the synthetic data generated by independent simulations. The presence or absence of connectivity in each direction is decided by whether or not an output value z ∈ [0, 1] exceeds a threshold θ . It is possible to reduce the number of FPs by shifting the threshold θ to a high level. But this operation may produce many FNs, making many existing connections be missed. To balance the false-positives and false-negatives, we considered maximizing the Matthews correlation coefficient (MCC) 44 , as has been done in our previous study 14 . The MCC is defined as where N TP , N TN , N FP , and N FN represent the numbers of true positive, true negative, false positive, and false negative connections, respectively.
We have obtained two coefficients for excitatory and inhibitory categories and taken the macro-average MCC that gives equal importance to these categories (Macro-average) 45 , MCC = (MCC E + MCC I )/2 as we have done in the previous study 14 . In computing the coefficient for the excitatory category MCC E , we classify connections as excitatory or other (unconnected and inhibitory); for the inhibitory category MCC I , we classify connections as inhibitory or other (unconnected and excitatory). Here we evaluate MCC E by considering only excitatory connections of reasonable strength (EPSP > 0.1 mV for the MAT simulation and > 1 mV for the HH simulation).
We have confirmed that the Matthews correlation coefficient exhibits a wide peak at about θ ∼ 0.5 (Fig. 10), and accordingly, we adopted θ = 0.5 as the threshold.

A large-scale simulation of a network of MAT neurons. To obtain a large number of spike trains that
have resulted under the influence of synaptic connections between neurons, we ran a numerical simulation of a network of 1000 model neurons interacting through fixed synapses. Of them, 800 excitatory neurons innervate to 12.5 % of other neurons with EPSPs that are log-normally distributed 14,46-49 , whereas 200 inhibitory neurons innervate randomly to 25% of other neurons with IPSPs that are normally distributed.
Neuron model. As for the spiking neuron model, we adopted the MAT model, which is superior to the Hodgkin-Huxley model in reproducing and predicting spike times of real biological neurons in response to fluctuating inputs 21,23 . In addition, its numerical simulation is stable and fast. The membrane potential of each neuron obeys a simple relaxation equation following the input signal: where g e , g i represents the excitatory conductance and the inhibitory conductance, respectively. Here RI bg represent the background noise. The conductance evolves with the where τ s,X is the synaptic time constant, X stands for e (excitatory) or i (inhibitory), t jk is the kth spike time of jth neuron, d j is a synaptic delay and G j is the synaptic weight from jth neuron. δ(t) is the Dirac delta function. Next, the adaptive threshold of each neuron θ(t) obeys the following equation: where t j is the jth spike time of a neuron, ω is the resting value of the threshold, τ k is the kth time constant, and α k is the weight of the kth component. The parameter values are summarized in Table 4.
Synaptic connections. We ran a simulation of a network consisting of 800 pyramidal neurons and 200 interneurons interconnected with a fixed strength. Each neuron receives 100 excitatory inputs randomly selected from 800 pyramidal neurons and 50 inhibitory inputs selected from 200 interneurons. The excitatory and inhibitory synaptic connections were sampled from respective distributions so that the resulting EPSPs and IPSPs are similar to the distributions adopted in our previous study 14 . In particular, the excitatory conductances {G E ij } were sampled independently from a log-normal distribution 46,47 . where µ = −5.543 and σ = 1.30 are the mean and SD of the natural logarithm of the conductances.
The inhibitory conductances {G I ij } were sampled from the normal distribution:   Background noise. Because our model network is smaller than real mammalian cortical networks, we added a background current to represent inputs from many neurons, as previously done by Destexhe et al. 11,50 .
The summed conductance RI bg represents random bombardments from a number of excitatory and inhibitory neurons. The dynamics of excitatory or inhibitory conductances can be approximated as a stationary fluctuating process represented as the Ornstein-Uhlenbeck process 51 , where g X stands for g e or g i , and ξ(t) is the white Gaussian noise satisfying �ξ(t)� = 0 and �ξ(t)ξ(s)� = δ ij δ(t − s). The real biological data has a wide variety of fluctuation, including non-trivial large variations with some characteristic timescales. For instance, the hippocampal neurons are subject to the theta oscillation of the frequency range of 3 − 10 (Hz) 52 . To reproduce such oscillations that are also observed in the cross-correlogram, we introduced slow oscillations into the background noise for excitatory neurons, as where ξ 1 (t) and ξ 2 (t) are the white Gaussian noise satisfying �ξ i (t)� = 0 and �ξ i (t)ξ j (s)� = δ ij δ(t − s).
Among N = 1000 neurons, we added such oscillating background signals to three subgroups of 100 neurons (80 excitatory and 20 inhibitory neurons), respectively with 7, 10, and 20 Hz. The phases of the oscillation δ were chosen randomly from the uniform distribution. Amplitudes of the oscillations were chosen randomly from uniform distribution in an interval [Ã/2, 3Ã/2] . The parameters for the background inputs are summarized in Table 5.
Numerical simulation. Simulation codes were written in C++ and parallelized with OpenMP framework. The time step was 0.1 ms. The neural activity was simulated up to 7200 s. Experimental data. Spike trains were recorded from the PF, IT, and V1 cortices of monkeys in three experimental laboratories using the Utah arrays. All the studies were carried out in compliance with the ARRIVE guidelines. Individual experimental settings are summarized as follows.
Prefrontal cortex (PF). The experiments were carried out on an adult male rhesus macaque Macaca mulatta (6.7 kg, age 4.5 y). The monkey had access to food 24 h a day and earned liquid through task performance on testing days. Experimental monkeys were socially pair housed. All experimental procedures were performed in accordance with the ILAR Guide for the Care and Use of Laboratory Animals and were approved by the Animal Care and Use Committee of the National Institute of Mental Health (U.S.A.). Procedures adhered to applicable United States federal and local laws, including the Animal Welfare Act (1990 revision) and applicable Regulations (PL89544; USDA 1985) and Public Health Service Policy (PHS2002). Eight 96-electrode arrays (Utah arrays, 10 × 10 arrangement, 400 μm pitch, 1.5 mm depth, Blackrock Microsystems, Salt Lake City, U.S.A.) were implanted on the prefrontal cortex following previously described surgical procedures 53 . Briefly, a single bone  www.nature.com/scientificreports/ flap was temporarily removed from the skull to expose the PFC, then the dura mater was cut open to insert the electrode arrays into the cortical parenchyma. Next, the dura mater was closed and the bone flap was placed back into place and attached with absorbable suture, thus protecting the brain and the implanted arrays. In parallel, a custom-designed connector holder, 3D-printed using biocompatible material, was implanted onto the posterior portion of the skull. Recordings were made using the Grapevine System (Ripple, Salt Lake City, USA). Two Neural Interface Processors (NIPs) made up the recording system, one NIP (384 channels each) was connected to the 4 multielectrode arrays of one hemisphere. Synchronizing behavioral codes from MonkeyLogic and eyetracking signals were split and sent to each NIP. The raw extracellular signal was high-pass filtered (1 kHz cutoff) and digitized (30 kHz) to acquire single-unit activity. Spikes were detected online and the waveforms (snippets) were stored using the Trellis package (Grapevine). Single units were manually sorted offline using custom Matlab scripts to define time-amplitude windows in combination with clustering methods based on PCA feature extraction. Further details about the experiment can be found elsewhere 54 . Briefly, the recordings were carried out while the animals were comfortably seated in front of a computer screen, performing left or right saccadic eye movements. Each trial started with the presentation of a fixation dot on the center of the screen and the monkeys were required to fixate. After a variable time (400-800 ms) had elapsed, the fixation dot was toggled off and a cue (white square, 2 • × 2 • side) was presented either to the left or right of the fixation dot. The monkeys had to make a saccade towards the cue and hold for 500 ms. 70% of the correctly performed trials were rewarded stochastically with a drop of juice (daily total 175-225 mL). Typically, monkeys performed > 1000 correct trials in a given recording session for recording time of 120-150 min.
Inferior temporal cortex (IT). The experiments were carried out on an adult male Japanese monkey (Macaca fuscata, 11 kg, age 13 y). The monkey had access to food 24 h a day and earned its liquid during and additionally after neural recording experiments on testing days. The monkey was housed in one of adjoining individual primate cages that allowed social interaction. All experimental procedures were approved by the Animal Care and Use Committee of the National Institute of Advanced Industrial Science and Technology (Japan) and were implemented in accordance with the "Guide for the Care and Use of Laboratory Animals" (eighth ed., National Research Council of the National Academies). Four 96 microelectrode arrays (Utah arrays, 10 × 10 layout, 400 μm pitch, 1.5 mm depth, Blackrock Microsystems, Salt Lake City, USA) were surgically implanted on the IT cortex of the left hemisphere. Three arrays were located in area TE and the remaining one in area TEO. Surgical procedures were roughly the same as having been described previously 53 , except that a bone flap that was temporarily removed from the skull was located over the IT cortex and that a CILUX chamber was implanted onto the anterior part of the skull protecting connectors of the arrays. Recordings of neural data and eye positions were done in a single session using Cerebus TM system (Blackrock Microsystems). The extracellular signal was bandpass filtered (250-7.5 k Hz) and digitized (30 kHz). Units were sorted online before the recording session for the extracellular signal of each electrode using a threshold and time-amplitude windows. Both the spike times and the waveforms (10 and 38 samples, preceding and after a threshold crossing, respectively) of the units were stored using Cerebus Central Suite (Blackrock Microsystems). Single units were refined offline by hand using the PCA projection of the spike waveforms in Offline sorter TM (Plexon Inc., Dallas, USA). The monkey seated in a primate chair, and the head was restrained with a head holding device so that the eyes were positioned 57 cm in front of a color monitor's display (GDM-F520, SONY, Japan). The display subtended a visual angle of 40 • × 30 • with a resolution of 800 × 600 pixels. A television series on animals (NHK's Darwin's Amazing Animals, Asahi Shimbun Publications Inc., Japan) was shown on the display throughout the online spike sorting and the recording session. The monkey's eye position was monitored using an infrared pupil-position monitoring system 55 and was not restricted.
The primary visual cortex (V1). The data set was obtained from Collaborative Research in Computational Neuroscience (CRCNS), pvc-11 56 by the courtesy of the authors of 57 . In this experiment, spontaneous activity was measured from the primary visual cortex while a monkey viewed a CRT monitor ( 1024 × 768 pixels, 100 Hz refresh) displaying a uniform gray screen (luminance of roughly 40 cd/m 2 ). Briefly, the animal was premedicated with atropine sulfate (0.05 mg/kg) and diazepam (Valium, 1.5 mg/kg) 30 min before inducing anesthesia with ketamine HCl (10.0 mg/kg). Anesthesia was maintained throughout the experiment by a continuous intravenous infusion of sufentanil citrate. To minimize eye movements, the animal was paralyzed with a continuous intravenous infusion of vecuronium bromide (0.1 mg/kg/h). Vital signs (EEG, ECG, blood pressure, end-tidal PCO2, temperature, and lung pressure) were monitored continuously. The pupils were dilated with topical atropine and the corneas protected with gas-permeable hard contact lenses. Supplementary lenses were used to bring the retinal image into focus by direct ophthalmoscopy and later adjusted the refraction further to optimize the response of recorded units. Experiments typically lasted 4-5 d. All experimental procedures complied with guidelines approved by the Albert Einstein College of Medicine of Yeshiva University and New York University Animal Welfare Committees. Spike sorting and analysis criteria: Waveform segments were sorted off-line with an automated sorting algorithm, which clustered similarly shaped waveforms using a competitive mixture decomposition method 58 . The output of this algorithm was refined by hand with custom time-amplitude window discrimination software (written in MATLAB; MathWorks) for each electrode, taking into account the waveform shape and interspike interval distribution. To quantify the quality of the recording, the signal-to-noise ratio (SNR) of each candidate unit was computed as the ratio of the average waveform amplitude to the SD of the waveform noise [59][60][61] . Candidates that fell below an SNR of 2.75 were discarded as multiunit recordings.