Tuning dissimilarity explains short distance decline of spontaneous spike correlation in macaque V1

Fast spike correlation is a signature of neural ensemble activity thought to underlie perception, cognition, and action. To relate spike correlation to tuning and other factors, we focused on spontaneous activity because it is the common 'baseline' across studies that test different stimuli, and because variations in correlation strength are much larger across cell pairs than across stimuli. Is the probability of spike correlation between two neurons a graded function of lateral cortical separation, independent of functional tuning (e.g. orientation preferences)? Although previous studies found a steep decline in fast spike correlation with horizontal cortical distance, we hypothesized that, at short distances, this decline is better explained by a decline in receptive field tuning similarity. Here we measured macaque V1 tuning via parametric stimuli and spike-triggered analysis, and we developed a generalized linear model (GLM) to examine how different combinations of factors predict spontaneous spike correlation. Spike correlation was predicted by multiple factors including color, spatiotemporal receptive field, spatial frequency, phase and orientation but not ocular dominance beyond layer 4. Including these factors in the model mostly eliminated the contribution of cortical distance to fast spike correlation (up to our recording limit of 1.4mm), in terms of both 'correlation probability' (the incidence of pairs that have significant fast spike correlation) and 'correlation strength' (each pair's likelihood of fast spike correlation). We suggest that, at short distances and non-input layers, V1 fast spike correlation is determined more by tuning similarity than by cortical distance or ocular dominance.


Introduction
Fast spike correlation ('synchrony' or 'coincident spiking') is thought to be a code or signature of spiking ensembles, enabling the brain to perform efficient computations relating to perception and behavior (Salinas & Sejnowski, 2001;Singer, 1999). Although the functional roles of different timescales of spike correlation is unclear, synchrony within a narrow time window, approximating the temporal integration window of downstream neurons acting as coincidence detectors, is considered separately from slower changes in correlated excitability (noise correlation, 'Rsc', (Cohen & Kohn, 2011)). To develop computational models of spiking ensembles, it is necessary to know which neurons fire coincidently (including before stimulus onset), and how much of the coincident firing is related to the neurons' tuning properties vs. due to cortical distance (Masquelier & Thorpe, 2007).
Although several reports have focused on the dynamic or context-dependent nature of spike correlations (Das & Gilbert, 1999;Gray et al., 1989;Hung, Ramsden, & Roe, 2007;Roelfsema et al., 1997;Stettler et al., 2002), here we focus on spike correlations during spontaneous activity because the variation in correlation strength across cell pairs is typically many times larger than the variation across stimulus conditions in the same pairs (Hung,  Ramsden, & Roe, 2011;Luczak, Bartho, & Harris, 2009), suggesting that the mechanism is mostly intrinsic (i.e. tied to the functional architecture) rather than stimulus-dependent (Ringach, 2009).
Although visual stimulation typically reduces neural variability (Churchland et al., 2010), correlations presented in both evoked and spontaneous activity are thought to share similar mechanisms (Jermakowicz et al., 2009). Understanding variability of spike correlations during spontaneous activity (i.e. whether this variability is random or systematic) is thus a necessary component to understanding signal processing in basic cortical circuits. Also, identifying the factors underlying spontaneous correlations may provide a fairer and more consistent baseline for comparing across studies that examine different stimuli and different mixtures of cells.
Although new methods combining slice physiology, in vivo calcium imaging, and/or electron microscopy have been developed that have enabled precise alignment between tuning, morphology, and circuitry in rodents and small animals (Bock et al., 2011;Ko et al., 2011;Lefort et al., 2009;Shepherd et al., 2005), such methods are extremely difficult to implement, especially in larger animals such as macaque monkeys which are more similar to humans. Analyses of spontaneous extracellular spiking may thus offer an enormous advantage in studying coordinated assembly activity.
Although many studies have suggested that coincident spiking declines with cortical distance, it is unclear whether this distance dependency is simply due to examining too few factors. It is well known that different tuning factors are related, and that the interaction of these factors may lead to a residual effect of cortical distance when only single factors are examined (e.g. a study may find that spike correlation depends on both orientation and distance, but the distance dependency may be due to an unexamined factor that is co-linear with distance). However, no study has sampled sufficient tuning properties (typically no more than 2 or 3, analyzed separately) and cell pairs to disentangle the effect of cortical distance from the effect of the overall decline in tuning similarity across combinations of tuning properties. Also, rather than sampling one site per penetration, it would be better to sample multiple sites per penetration (multiple pairs of neurons with equal horizontal separation and the same topography) to disentangle this relationship.
Here, we asked what are the relative contributions of different tuning properties to spike correlation, and whether horizontal cortical distance has a separate contribution, beyond that already predicted by tuning dissimilarity. The standard hypothesis is that spike correlation depends on horizontal cortical distance (Das & Gilbert, 1999;Gray et al., 1989;Hata et al., 1991;Hung, Ramsden, & Roe, 2007;Maldonado, Friedman-Hill, & Gray, 2000;Smith & Kohn, 2008;Toyama, Kimura, & Tanaka, 1981), in addition to tuning similarity (e.g. for orientation, ocular dominance, chromatic preference, and spatiotemporal receptive field similarity (Das & Gilbert, 1999;DeAngelis et al., 1999;Engel et al., 1990;Hata et al., 1991;Nowak et al., 1995;Schwarz & Bolz, 1991;Ts'o & Gilbert, 1988;Ts'o, Gilbert, & Wiesel, 1986)). However, whether if and to what extent spike correlations actually depend on cortical distance is unclear, because not all factors were significant (Chiu & Weliky, 2002), including horizontal cortical distance (Samonds et al., 2006;Schwarz & Bolz, 1991), and horizontal interactions can be found at up to 4-7 mm (Engel et al., 1990;Smith & Kohn, 2008), and even across hemispheres (Bosking et al., 2000;Engel et al., 1991;Nowak et al., 1999). We suggest an alternative hypothesis, that the decline of fast spike correlation with cortical distance, at least for short distances (<1.4 mm, about 1-2 hypercolumns), can be explained by the decline in tuning similarity with distance (this possibility was also mentioned in (Ts'o, Gilbert, & Wiesel, 1986)). If so, it should be possible to use GLM to precisely quantify weights (beta coefficients) for both spike correlation probability and correlation strength, and to determine whether cortical distance is a significant predictor beyond that already predicted by tuning similarity.

Animal preparation and surgery
We recorded from two 4-5 kg Formosan macaque monkeys (Macaca cyclopis). M. cyclopis is a member of the group M. mulatta along with M. fuscata, and is paraphyletic to M. nemestrina and M. fascicularis based on mitochondrial DNA sequences (Li & Zhang, 2005). All experimental procedures were performed in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and were approved by the Institutional Animal Care and Use Committee of National Yang-Ming University.
Anesthesia was induced with ketamine (10 mg/kg IM). Monkeys were artificially respired and continuously monitored for EEG, EKG, body temperature, expired CO 2 , and pO 2 . Light anesthesia was maintained with sodium thiopental (Pentothal 2 mg/kg/h IV) and the muscle relaxant rocuronium bromide (Esmeron 1.2 mg/kg/h IV), and anesthetic depth was maintained via custom software that continuously measured delta vs. gamma EEG power. Pupils were dilated with atropine sulfate. The center of gaze was estimated via reverse ophthalmoscopy of the optic discs, and the eyes were focused via contact lenses and converged upon the monitor at 57 cm distance.
At the start of each session, the eyes were converged via rotating wedge prisms (Thorlabs RSP1). Based on the alignment of small $0.2 deg receptive fields from the two eyes (three neurons per monkey), we estimate the precision of alignment to be less than 0.1 deg. Except for ocular dominance (OD) measurements, all other blocks were presented monocularly to avoid the possibility of phase mis-alignment of receptive fields from the two eyes (phase mis-alignment alone could elevate OD as a factor). The longest block, spatiotemporal receptive field (STRF) mapping, took 25 min and yielded RFs as small as 0.44°wide with 0.22°wide subfields, indicating that the eyes were stable throughout the recording. The ''orientation/SF/phase'' recording block took 15 min and also yielded reliable phase preference even at 2 cyc/deg ( Fig. 2E and F).

Electrophysiology
We inserted 64-site multi-electrode arrays (A8 Â 8-5 mm200-200-413, 8 penetrations ('shanks'), 8 sites per penetration, spanning 1.4 Â 1.4 mm horizontally and in depth, 200 lm spacing, Neuronexus Technologies, Inc.) normal to the cortical surface, 14 mm anterior of the occipital ridge and 10 mm lateral of midline (approximately 3-4 deg eccentricity). The width of the array was thus sufficient to span two complete cycles of ocular dominance hypercolumns (as measured in a third monkey by aligning the array across OD columns). Cortical depth was assessed by DiI and cytochrome oxidase staining ( Fig. S1A and B), current source density analysis (Fig. S1C) and by the temporal frequency limit outside layer 4 ( Fig. S1D-H). Spikes (400-5000 Hz) and local field potentials (LFPs, 1-300 Hz) were filtered (48 dB/octave) and continuously digitized at 24.4 kHz (RZ2, Tucker-Davis Technologies, Inc.). Single units were isolated offline via super-paramagnetic clustering (Quiroga, Nadasdy, & Ben-Shaul, 2004). To avoid possible errors from the unsupervised spike sorting algorithm, we rejected and manually resorted all spike clusters ('units') if over 5% of interspike intervals were <2.5 ms. Manual sorting was done by adjusting the temperature of the annealing in the super-paramagnetic clustering algorithm. At low temperature, all spikes are assigned to the same cluster, whereas at high temperature, each spike forms a single cluster. We chose an 'optimal' temperature by gradually increasing the temperature until less than 5% of interspike intervals were <2.5 ms. This criterion is considered 'good' in extracellular recordings (e.g. (Tolias et al., 2007) had 8% falsely assigned spikes). We matched the same neuron across different recording blocks by choosing the waveforms that had the highest correlation. Because typically each channel only had 1 or 2 significant waveforms and the match correlations were very high (r > 0.97, vs. $0.62 for nonmatch), the matching was not difficult.
Multiple detection of the same unit across different channels was tested via cross-correlation analysis with 1 ms histogram bins (Hung, Ramsden, & Roe, 2007). Of 3503 pairs tested, 101 (2.8%) were excluded as possible cases of multiple detection because they had narrow 1 ms peaks, defined by B 0 > (B 1 + B À1 ) and B 0 > 0.05, where B 0 is the height of the bin (in units of coincidences per spike pair) at 0, À1, or 1 ms (any of the 3) and B 1 , and B À1 are the heights of the neighboring bins.

Visual stimulation
Stimuli were generated in Matlab (Mathworks Inc.) and Psychophysics Toolbox 3 (Brainard, 1997) and displayed on a 21-in. calibrated CRT monitor (ViewSonic pf227). Phosphor emission spectra were measured via spectrometer, and stimuli were gamma corrected. Following hand mapping of receptive fields (RFs), all stimuli were centered on the cluster of RFs. The most distant receptive fields were less than 1°apart. All stimuli were presented at an average luminance of $27.8 cd/m 2 , about half the maximum luminance of our monitor (Commission Internationale de l'Eclairage 1931; x, 0.3106; y, 0.3261). Ocular dominance was measured via full contrast chromatic drifting squarewave gratings (red, green, blue, and yellow interleaved with black; 4 Hz, 2 color/black cycles per deg, i.e. 0.25 deg per stripe, 4 orientations). Random monocular presentation to left or right eye was controlled by eye shutters, with shutters closed between presentations. Other stimulus dimensions were tested via monocular presentation to the contralateral eye.
Spatiotemporal receptive fields were determined by spike-triggered analysis. Spike-triggered analysis was based on spike-triggered average (STA) and spike-triggered covariance (STC) analysis of random checkerboard patterns preceding each spike. Stimuli were binary checkerboard m-sequence patterns extending 7.04°v ertically and horizontally, a 32 Â 32 grid of 0.22°Â 0.22°patches, updated at 42.5 Hz. We placed light anchors (4 Â 4 deg, full luminance) at the four corners of the screen to stabilize brightness (Gilchrist et al., 1999). Chromatic and luminance tuning were measured via flashed homogenous color patches (5 Hz, 4 Â 4 deg square) (Wachtler, Sejnowski, & Albright, 2003). The colors consisted of 8 hues with equal psychometric distance (based on the German Standard Color Chart) and 4 luminance levels (39.3, 33, 27.8, 23.4 cd/m 2 ). -stimulus period and ranges from 1 (ipsilateral) to À1 (contralateral). The distribution was broad at each cortical depth (from 0.2 mm down to 1.6 mm, boxes show quartiles and extremes), with a bias towards contralateral units at 1.0-1.2 mm, consistent with monocular input to layer 4 and array alignment along ocular dominance columns. (F) Distribution of OD similarities (N = 3503 pairs across two monkeys). OD similarity between each cell pair was defined as 1 À abs(OD_index_1 À OD_index_2). Arrows indicate similarities for these example pairs.

Selection criterion for neurons in GLM
The neurons were located across all cortical layers and most recording sites (at least 41 out of 64 sites per monkey). We measured neurons with at least 500 spikes during 15 min of spontaneous activity (i.e. >0.56 Hz). This minimum spike count threshold was met by 107/117 units in M1 and 62/78 units in M2 (169 total units), yielding 5671 pairs in M1 and 1891 pairs in M2. This set was further reduced to 118 units (68 in M1, 50 in M2) that had significant receptive fields based on spike-triggered analysis (see spatiotemporal receptive field (STRF) similarity), yielding 3503 cell pairs (2278 in M1, 1225 in M2) in the final analysis.

Identification of cortical layers
Identification of cortical layers in V1 was based on several methods. Cytochrome oxidase (CO) staining in M3 indicated the location of color-sensitive blobs that mark output layers (Livingstone & Hubel, 1984;Lu & Roe, 2008). Current source density (CSD) analysis in M1 was calculated from the trial averaged LFP evoked with static gratings averaged across all orientations. CSD was approximated by a 3-point formula to the second spatial derivative of the LFP: where h(n) is the LFP signal on the nth electrode contact, h(n À 1) and h(n + 1) is the LFP signal of the immediately adjacent contact above and below the nth contact on the same shank, and D is the spacing between contacts (0.2 mm). Layer 4 was determined by the reversal of the CSD compared with superficial or deep layers (Kajikawa & Schroeder, 2011;Mitzdorf, 1985).
In addition, we tested gratings at several temporal frequencies in V1 (1.7, 3.5, 7, 14, 28 Hz) and found a difference in the high temporal frequency cutoff between input and output layers, consistent with a previous report of a decrease of high-cutoff temporal frequency of about 20 Hz between LGN and V1 output layers (Hawken, Shapley, & Grosof, 1996).

Spike correlation
Spike correlation between cells was measured via crosscorrelation analysis (Ts'o, Gilbert, & Wiesel, 1986) from continuous periods of spontaneous activity (15 min with the screen off). To avoid the possibility of plasticity from stimulus presentation (Li & DiCarlo, 2010;Yao & Dan, 2001), spike correlations during spontaneous activity were measured in a separate recording block prior to measuring tuning properties. To avoid bias from spike collisions (missing spike correlation on the same channel because they occur too closely or are overlapping in time), we only analyzed pairs recorded from different channels, not pairs sorted from the same channel.
The 'raw' cross correlation histogram (CCH) was defined as where S i and S j are binary values representing neuron 1's ith spike and neuron 2's jth spike at a given temporal offset s (s's are binned in the CCH at 1 ms per bin), and their product indicates whether an instance of a specific timing difference occurred (i.e. where neuron 1's spike preceded neuron 2's spike by s). s is limited by the size of the convolution window (from À150.5 to 150.5 ms at 1 ms per bin, total 301 bins). The total is normalized by the number of spikes in the two spike trains (N 1 and N 2 ) so that recording longer or shorter durations does not alter the CCH. Although previous studies corrected spike correlations for stimulus-locked responses, they only recently began controlling for slow correlations in spike timing (e.g. slow changes in excitability that are driven by modulatory neurotransmitters, traveling waves, or large-scale coupling between distant areas) that are not stimuluslocked and that can distort the measurement of spike correlation (Brody, 1998;Sato, Nauhaus, & Carandini, 2012;Smith & Kohn, 2008). These limitations make it difficult to determine how different factors relate to fast spike correlation, necessary to construct computational models of ensemble dynamics (Carandini & Ringach, 1997;Dehaene & Changeux, 2005;Goldberg, Rokni, & Sompolinsky, 2004;Grossberg & Williamson, 2001;Li, 1998;Masquelier & Thorpe, 2007;McLaughlin et al., 2000;Riesenhuber & Poggio, 2000).
We corrected the CCH for firing rate correlations slower than 50 ms by using the jitter method (Smith & Kohn, 2008). This jitter window size was at the narrow end of those tested by Smith and Kohn. This correction is more conservative than standard shift or shuffle predictors, because it is based on spontaneous activity and because it removes covariations in firing rate slower than 50 ms (however, it cannot remove all sources of common input).
We measured the significance of the spike correlation based on the [À50.5:50.5] ms interval. We permuted the spike train 2020 times via the jitter method and defined the significance threshold ('95%ile jitter predictor') as the maximum (per bin) of the 2020 permutations. We used 2020 permutations because there are 101 bins in the [À50.5:50.5] ms interval, and so the likelihood that any bin in the peak exceeds the significance threshold is 1/2020 bin À1 Â 101 bins (Bonferroni correction) = 101/2020, i.e. p = 0.05. We defined correlation significance as a binary value indicating whether the CCH passed the threshold.
We defined 'correlation probability' as the incidence of pairs that had significant correlation (i.e. the total number of significant pairs, normalized by the number of tested pairs). We defined 'correlation strength' as the height of the jitter-corrected correlogram (i.e. CCH minus the 50%ile of the jitter predictor), also in the [À50.5:50.5] ms interval, i.e. the strength of the spike correlation that is faster than 50 ms.

Ocular dominance similarity
Ocular dominance was measured in response to chromatic drifting squarewave gratings presented for 2 s 'On' and 2 s 'Off' (see Section 2.3). For each cell, we calculated the 'evoked' (baseline-subtracted) responses in the [0.05:2.05] s interval after stimulus onset, averaged across 10 repetitions. The baseline was the average response in the [À150:À50] ms interval across all conditions, i.e. one number per cell, not stimulus-specific. The ocular dominance index was defined as where rI and rC are the evoked responses to the ipsilateral and contralateral eyes. We then defined ocular dominance similarity (ODS) as follows: where A and B are the two cells. ODS ranges from 1 to À1. An ODS of 1 indicates that both cells have exactly the same eye preference, whereas an ODS of À1 indicates that the cells have opposite eye preferences (and that both are monocular).

Spatial frequency, phase, and orientation similarity
We measured each cell's evoked (baseline-subtracted) responses (mean of [40:70] ms) to phase, spatial frequency, and orientation, averaged across 20 repetitions of flashed achromatic sinusoidal gratings (94 ms On, 106 ms Off presentation, see Section 2.3 and Fig. 2A, C, and E) and converted this 8 Â 8 Â 4 array (8 spatial frequencies Â 8 phases Â 4 orientations) to a 1 Â 256 array. The baseline was the [0:30] ms interval averaged across all conditions (one number per cell, not stimulus-specific). We then defined R signal of each pair, i.e. their tuning similarity in terms of spatial frequency, phase, and orientation, as their Pearson correlation.

Spatiotemporal receptive field (STRF) similarity
We measured each cell's spatiotemporal receptive field (STRF) using spike-triggered analysis. Its advantage is that it makes fewer assumptions about the cell's preferred stimulus and response patterns. We presented random binary checkerboard patterns (generated via m-sequence, (Reid, Victor, & Shapley, 1997)). The STRF was measured as spike-triggered average (STA), which captures the average pixel patterns preceding each spike, and as spike-triggered covariance (STC), which captures pixel covariation patterns preceding each spike.
Spike-triggered analysis estimates the low dimensional linear subspace of the full stimulus space (Chichilnisky, 2001;de Ruyter van Steveninck & Bialek, 1988;Paninski, 2003;Schwartz et al., 2006) and is similar to the Wiener/Volterra approach to determine the firing rate function of a neuron (Korenberg, Sakai, & Naka, 1989). Each spike provides a reference point for its preceding stimulus frames. We averaged all the spike-triggered stimulus frames in a defined temporal window preceding the spike across the entire stimulus space (spike triggered average, STA) as follows, where N is the total number of spikes recorded, and S i is the stimulus preceding the ith spike. We used a spike-triggered analysis window of 141 ms (23.5 ms per frame Â 6 frames). In previous studies, STA was used to map specific cone inputs to LGN receptive fields (Reid & Shapley, 2002) and receptive fields of simple cells in the primary visual cortex (DeAngelis, Ohzawa, & Freeman, 1993;Jones, Stepnoski, & Palmer, 1987). We used spike-triggered covariance (STC) to capture additional filters (low-dimensional subspaces), an approach that has been used in many systems to characterize multidimensional models via nonlinear combination rules (Rust et al., 2005;Touryan, Lau, & Dan, 2002). The spike-triggered covariance matrix, which embodies the multidimensional covariance structure of a neuron's receptive field, is computed as follows: By applying eigenvalue decomposition, we identified significant eigenvectors from the STC matrix. The first eigenvector is the direction in stimulus space along which the spike-triggered stimuli have the greatest variance, and subsequent eigenvectors are orthogonal directions in stimulus space explaining successively lower variance. Typically, only the first one or two eigenvectors were significant (exceeding the 95%ile of the explained variance of the shuffled distribution), based on shuffling the stimulus sequence (we ignored STCs higher than 2, see Selection criterion for neurons in GLM). The STRF thus may contain 0 or 1 significant STAs and 0 to 2 STCs.
We measured the similarity between STRFs in several ways. A simple correlation coefficient was computed between the STAs of each pair of cells (DeAngelis et al., 1999). For each cell's STA, we concatenated the three consecutive frames (frame 2-4, i.e. (23.5:94] ms preceding the spike) that had the greatest spatial variance.
To compute the correlation between the STA of one cell vs. the STC of another cell, we applied two methods, one based on RF energy and the other based on the STC kernel itself or its negative (for cells with only one significant STC) or Procrustes transformation of the STC kernel (for cells with two significant STCs). Because the RF energy and Procrustes transformation methods yielded similar results, we report only the RF energy results.
The RF energy method was based on the correlation between the energy of receptive fields (receptive field overlap). Receptive field energy is defined as follows: where RF could be any significant STA or STC kernel.
We tried several variations of computing a summary score of STRF correlation, e.g. by taking the maximum Pearson correlation value across these methods (i.e. Max of STA vs. STA, STA vs. STC, and STC vs. STC comparisons). However, in reporting GLM results, the functional similarity of STRF was based on only the STA-STA similarity measure, because it explained better than STA-STC and STC-STC combinations, including RF energy (RF overlap, see results).

Full field color similarity
We flashed homogenous color patches at 5 Hz (94 ms ON, 106 ms OFF, (Wachtler, Sejnowski, & Albright, 2003)). For each cell, we measured color responses across 32 conditions (8 hues Â 4 luminances) and converted these evoked responses (mean of [40:130] ms, averaged across 20 repetitions, minus baseline in the [0:30) ms interval, one baseline value per cell) to a 1 Â 32 array. We then defined the color tuning similarity of each cell pair as the Pearson correlation of their arrays.

Generalized linear models
To quantify the relationship between spike correlation (correlation probability or correlation strength) and tuning similarity, we developed a generalized linear model (GLM) with the tuning similarity measurements and cortical distance as factors.
For correlation probability, we used a logistic regression form of the GLM, because correlation significance is a binary (not continuous) measure of the presence (1) or absence (0) of spike correlation faster than 50 ms. We began with a column vector Z, consisting of binary numbers Z i , i = 1,2,. . ., N, indicating whether a particular cell pair i had a significant correlogram peak. This was matched to a design matrix X, composed of N rows and K + 1 columns, where K is the number of factors in the model. For each row of the design matrix, the first element x i0 = 1. Logistic GLM was defined as follows: where f(x) is computed for each pair as a continuous variable, using beta coefficients b k learned from the training population. p i = P(Z i = 1) estimates the probability that a particular pair i has significant correlation. x ik is cell pair i's tuning similarity for factor k, b k is the coefficient for that factor, with K total factors and N total pairs. b 0 is the offset for each monkey. The large difference in b 0 across monkeys appears to be related to a difference in average firing rate (9.9 sp/s for M1, 4.3 sp/s for M2) and does not appear to affect the beta coefficients. The intuition behind this regression is to use a linear combination f(x) of multiple tuning similarity measures (x) to predict whether a pair's spike correlation is significant (Z = 1) or not (Z = 0). Fig. 9A and C shows the relationship between the data and the fitted model. The pairs are sorted along the abscissa according to f(x). Although Eq. (1) (solid lines) was fitted to the binary (peak) significance data, because the binary data are difficult to see, we plot the binary significance data of nearby pairs by applying a smoothing window of 15 pairs (dotted lines) to the binary data.
For correlation strength (the jitter-subtracted correlogram peak height of significantly correlated pairs), we used a logarithmic regression form of the GLM because the peak height distribution was normally distributed after logarithmic transformation (as verified via quantile-quantile plots and Komolgorov-Smirnov testing). The GLM for correlation strength was as follows: where Y i is the correlation strength for one cell pair, x ik is cell pair i's functional tuning similarity for factor k, b k is the coefficient for that factor, with K total factors and N total pairs. b 0 is the offset for each monkey. For computing R 2 during cross-validation (Eq. (6)), p i is de-

Model selection
We screened factors to include in the combined model by identifying factors that were significant at p < 0.2 with single-factor models. Screened factors were then included in our final (reduced) model. We then sequentially reduced the number of factors in the model, by removing the factor with the highest p-value, until only significant (p < 0.05) factors remained. Finally, because some factors were co-linear, we used the log likelihood ratio test to verify that the reduced model was not significantly less likely than the full model.
The standard equation for computing likelihood (L) for logistic regression is as follows: where y i is the number of pairs with significant correlogram peaks in a given population, and n i is the number of total pairs in that population. Because we only have a single population, n i = 1 and y i = Z i is a binary value (either 0 or 1, indicating whether a pair has a significant correlogram peak). The equation simplifies to the following: The likelihood (L) for logarithmic regression is as follows: 2.4.10. Evaluation of model performance To evaluate the GLM performance, we calculated explained variance (average R 2 ± SEM Â 100%) based on repeated cross-validations (N = 200) within each animal. We based our conclusions on cross-validation across cells because it is more conservative, but we also show the cross-validation across cell-pairs for completeness. For 'cell cross-validation', model coefficients were first estimated with the pairs formed by half of the cells and then tested on pairs among the remaining cells (i.e. no re-use of the same cells). For 'pair cross-validation', model coefficients were first estimated with 70% of total pairs and then tested on the remaining 30%.
More specifically, we estimated beta coefficients (b k ) and offset (b 0 ) using a subset of data (''training''). Based on those weights and offsets (learning the hyperplane), we evaluated the model performance on the remainder of the data (''testing''). The testing was based on Eq. (1) for binary logistic regression and Eq. (2) for logarithmic regression. The goodness of fit (GOF) of the logistic model was assessed with Hosmer-Lemeshow test (Hosmer et al., 1997).
Model performance was assessed via R 2 , which indicates the percent of variance in the data that is explained by the model. However, because R 2 is not frequently used in logistic regression, we adopted the definition of pseudo R 2 in Efron (1978), which has a similar formula as in linear regression. Both R 2 and pseudo R 2 were computed as follows: where Y i are the dependent variables in the regression model, p i are the predicted values based on Eqs. (1) and (2), N is the total number of pairs in the model testing. Because pseudo R 2 is not a standard measure of model performance of logistic regression, we also computed the ROC (receiver operating characteristic) curve. The curve is based on the true positive vs. false positive distribution between actual pairwise correlation significance (Y) and model predicted probability (p).
Model performance was then assessed as AUC (area under ROC curve). To determine the significance of the AUC vs. 0.5, we carried out a z-test where the z-score is the AUC minus expected chance 0.5, divided by the standard deviation of the AUC (across 20 resamples of cell or pair cross-validation).

Vertical (within-penetration) and horizontal (within-row) shuffling
To test whether the model results were due to a sampling bias due to the relationship between the array and the local topography (i.e. simply due to vertical (columnar) organization), we applied within-penetration shuffling, computed 20 cross-validations of each, then applied a t-test to the two distributions. 'Withinpenetration shuffling' was done by randomly replacing, for each pair, its correlation significance value (0 or 1) and correlation strength with that of a different pair of neurons. Each replacement neuron was chosen from the same or adjacent electrode shank (horizontal separation 6 0.2 mm) as its original and had average firing rate within ±5 sp/s of the original. As a separate control, we also shuffled 'within-row' (vertical separation 6 0.2 mm), which should further dilute the influence of topographic organization and result in AUC closer to chance (0.5).
These results are reported as p-values in Results: Correlation probability and correlation strength.

Results
We measured the tuning and spontaneous spike correlation of 118 V1 single units in two monkeys (68 in M1, 50 in M2), after selecting for cells that were active at >0.56 Hz and that had significant receptive fields based on spike-triggered analysis (see Sections 2 and 4). This yielded 3503 cell pairs (2278 in M1, 1225 in M2) in the final analysis. Spike correlation strengths were highly correlated between visually evoked and spontaneous conditions (M1: r = 0.67, p < 10 À10 ; M2: r = 0.54, p < 10 À10 , Pearson correlation of log spike correlation strengths, Fig. S2), indicating that variations in spike correlations are much larger across pairs than across conditions.

Ocular dominance
We presented full contrast color squarewave drifting gratings on separate trials to each eye. Cells responded vigorously and showed clear periodic response modulation. Three example cells ( Fig. 1A-C) show typical responses, preferring either one eye ( Fig. 1A and B) or both eyes (Fig. 1C). Grating responses appeared half-wave rectified ( Fig. 1A and C) or full-wave rectified (Fig. 1B). We quantified ocular dominance preference via a contrast measure (OD index) ranging from À1 (contralateral eye) to 1 (ipsilateral). OD indices were widely distributed and varied with cortical depth (Fig. 1D: M1; Fig. 1E: M2). OD was biased contralaterally in both monkeys, consistent with the alignment of our electrode arrays anterior-posteriorly along OD columns (guided by optical imaging in M1) to reduce the association between ocular dominance similarity and distance (we confirmed these results by aligning the array across OD columns in a third monkey, see below). The depth of maximum contralateral response was $1-1.2 mm, consistent with monocular input to layer 4 and with DiI, cytochrome oxidase staining, current source density analysis, and temporal frequency analysis (Fig. S1, see Section 2). OD similarity (ODS, see Methods) was skewed toward 1 (both cells had similar OD) but included a broad range of similarities across the two monkeys (Fig. 1F).

Spatial frequency/phase/orientation
We characterized the cells' preferred spatial frequency, phase, and orientation by flashing static Gabor gratings (see Section 2). Most cells showed transient and robust responses. Response patterns ranged between two types. Cells that exhibited simple celllike behavior were selective for different phases at different spatial frequencies, typically showing diagonal bands in the phase/frequency/orientation plot ( Fig. 2A, B, E and F). Cells that exhibited complex cell-like behavior were selective for a particular orientation and spatial frequency regardless of phase ( Fig. 2C and D). Fig. 2G shows the unimodal distribution of pairwise similarities (measured as Pearson correlation between baseline-subtracted response patterns; arrow indicates similarity for this pair).

Spatiotemporal receptive field (STRF)
We mapped each cell's STRF via spike-triggered averaging (STA, 'reverse correlation') and spike-triggered covariance (STC) analysis (Schwartz et al., 2006). Fig. 3A, C, and E shows STRFs of three example units (different cells than in Fig. 2). In all examples, the STRF pattern was strongest at $70.5 ms before the spike (stimulus frames are 23.5 ms apart, i.e. 2 video frames at 85 Hz), beginning as early as 47 ms.
The example cell in Fig. 4B and C was narrowly tuned for red across a range of luminances, whereas the example cell in Fig and E was broadly tuned to hue. Fig. 4F shows the distribution of hue preferences, based on cosine fitting of hue responses at each cell's preferred luminance. This distribution is consistent with previous reports (Wachtler, Sejnowski, & Albright, 2003), including a bias toward darker colors (Fig. 4G) (Yeh, Xing, & Shapley, 2009). Both chromatic and luminance histograms were similar across the two monkeys. Overall, 36/118 cells (31%) were tuned to hue, 38/118 cells (32%) to luminance, and 7/118 cells (6%) to the interaction of hue and luminance, ANOVA). As with other tuning properties, the distribution of color tuning similarities, based on Pearson correlation of hue and luminance responses, was unimodal and broad (Fig. 4H).

Horizontal distance
Horizontal cortical distance between cell pairs ranged from 0 mm (same penetration, different electrode contacts) to 1.4 mm (first and last penetrations) in 0.2 mm steps. On average, tuning similarity declined with distance for all functional properties except ocular dominance (OD), which did not vary with distance ( Fig. 5) because the arrays were aligned along ocular dominance columns. Because each distance had a wide variation in OD preference, we were still able to assess the effect of OD similarity via the GLM.

Conversion to tuning similarity indices for the GLM
Ocular dominance similarity was calculated as 1 minus the absolute difference of the OD indices (0-2) and so ranged from 1 (similar) to À1 (dissimilar) (Fig. 1F). Other functional similarities were measured as the Pearson correlation of the tuning properties, thus also ranging from 1 (similar) to À1 (dissimilar). Overall, these tuning similarity indices were unimodally distributed. A few were skewed due to the placement of the electrode array along OD columns (Fig. 1F) or due to the quantification of similarity (Figs. 2G and 3G).

Characterization of spike correlation
We used cross-correlation analysis to measure the precise spike timing differences of pairs of spike trains (Fig. 6A, black solid line: Raw cross-correlation histogram, 'CCH'). Recent reports have highlighted differences between fast and slow correlations (Benucci, Frazor, & Carandini, 2007;Kenet et al., 2003;Smith & Kohn, 2008;Xu et al., 2007). To focus on fast correlations that have been linked to both tuning similarity and cortical distance and that are more likely to trigger coincidence detection by downstream neurons, we used jitter analysis to estimate and remove correlations slower than 50 ms (Smith & Kohn, 2008). We report these jittercorrected results as correlation probability ('CP') and correlation strength ('CS'). CP is the incidence of pairs with significant spike correlation (see Section 2: Spike correlation). CS is measured for significant pairs only and is based on correlations faster than 50 ms, i.e. the height of the raw cross-correlation peak minus the 50%ile of the jitter predictor (see Section 2). Correlations slower than 50 ms are postulated to be due to non-specific factors like traveling waves or top-down feedback.
Example correlograms and rasters (black and gray dots, Fig. 6B-G) indicate that the sharp peaks were not simply due to artifacts of slowly covarying firing rate, e.g. due to slow covariations in excitability (Fig. 6C, E, and G). Although correlation significance and correlogram peak height are not entirely separate (weaker correlogram peak heights are also more likely to be non-significant), the transition was smooth across a wide range of peak heights (Fig. 7).
Decline in correlation probability with horizontal cortical distance. Consistent with previously reports, raw CCH peak height declined with horizontal cortical distance ( Fig. 6H and I). This decline in raw spike correlations contains a mixture of both fast and slow correlations. Correlation probability and correlation strength focus on the fast correlations that are postulated to be more focal and linked to tuning and possibly distance (Benucci, Frazor, & Carandini, 2007;Kenet et al., 2003;Smith & Kohn, 2008;Xu et al., 2007).
Correlation probability declined with distance in both monkeys ( Fig. 6J and K, for M1 and M2 resp., both R = À0.9, p < 0.001 based on mean correlation probability), whereas correlation strength weakly declined with distance in M1 (p = 0.02) but not M2 ( Fig. 6L and M). This weaker decline for correlation strength became non-significant upon cross-validation, and the decline for correlation probability became non-significant when receptive field tuning dissimilarity was included in the GLM (see below and Section 4).
The magnitudes of correlation probability differed across the two monkeys. We speculate that this is due to a difference in their sensitivity to anesthesia rather than a difference in sampling, because it coincided with a difference in the mean firing rates (9.9 ± 0.1 sp/s in M1, 4.3 ± 0.1 sp/s in M2). The higher firing rate of M1 produces higher jitter predictors, which elevates the significance threshold for calculating correlation probability. Despite these differences, the model beta coefficients and explained variances were similar across the two monkeys (see below).

Dependence of spike correlations on ocular dominance, color sensitivity, and cortical depth
Because the main analysis is based on measures of similarity, it does not specifically analyze interactions such as those between binocular vs. monocular cells, color sensitive vs. color insensitive cells, or interactions that are specific to layer 4. Fig. 8 shows how the number of significant pairs, correlation probability (CP), and correlation strength (CS) of significant pairs vary across a wide range of tuning interactions. As expected from the higher number of contralateral OD cells sampled, many significant pairs are contra-contra, with slightly more contra-binoc pairs than binoc-binoc pairs (Fig. 8A). CP corrects for this sampling bias by normalizing for the number of pairs tested for each type of interaction. CS examines significant pairs only. Neither CP nor CS varied strongly across the range of OD interactions tested. However, all three measures are clearly dominated by the contralateral OD input to layer 4, showing both contra-contra and contra-binoc interactions (Fig. 8B, G and L). We confirmed these results for OD similarity in a third monkey in which the array was placed across OD columns (Fig. S3).
The prevalence of color-insensitive ('gray') cells explains the higher number of significant pairs for gray-gray interactions (Fig. 8C). CP is higher for color-color pairs (Fig. 8H), and CS is similar across the range of color-gray interactions (Fig. 8M). These patterns are similar in layer 4 (Fig. 8D, I, and N). All three measures of spike correlation are higher for pairs of neurons at similar depths, particularly near layer 4 (Fig. 8E, J, and O).

Generalized linear model
We used a generalized linear model (GLM) to relate tuning similarity to spike correlation, using the pairwise similarity of each tuning property and horizontal distance as regressors (X) to explain spike correlation (Y or Z). We generated two types of GLMs, a logistic model for correlation significance (Z = 0 or 1) and a logarithmic model for peak height (Y) (see Section 2). In both cases, we primarily report results based on cell-wise cross-validation. Compared to cell-wise cross-validation, pair-wise cross-validation (i.e. training on some pairs, testing on other pairs, but allowing for re-use of the same cells) only slightly increased performance.  Fig. 9A and C shows the smoothed and predicted correlation probabilities based on logistic regression for the two monkeys. The GLM predicted about 13.6 ± 0.5% of the variance in correlation probability in M1 and 8.6 ± 0.3% in M2 (Efron's pseudo R 2 , see Section 2). Although correlation probability was higher in M2 than in M1, the goodness of fit (GOF) was similar across the two monkeys (v 2 11.38, p = 0.18 for M1; v 2 13.34, p = 0.10 for M2; non-significance indicates good fit, Hosmer-Lemeshow test for logistic regression (Hosmer et al., 1997).

Correlation probability
Because explained variance is not commonly used for logistic regression, we also measured performance based on signal detection theory. We plotted the receiver operating characteristic (ROC) curve for each monkey based on the binary correlation significance and predicted correlation probability (Figs. 9B and D). The area under the curve (AUC) was significantly higher than chance (M1: 0.69 and 0.70 for cell-and pair-wise cross-validation, respectively; M2: 0.63 and 0.65, resp.; p < 10 À5 in all cases, z-test vs. chance 0.5). It was also significantly higher than expected from columnar or layer-specific organization, based on vertical (withinpenetration) and horizontal (within-row) shuffling (M1: 0.57 and 0.53 for vertical and horizontal, resp.; M2: 0.54 and 0.52, resp.; p < 10 À10 in all cases, based on n = 20 cell-wise cross-validations, two-tailed t-test).
The model performance was not simply due to vertical (columnar) organization (i.e. not because neurons along the same penetration are similarly tuned or because of biased sampling of our array). Within-penetration shuffling (see Section 2) reduced the prediction to near chance (total explained variance = 0.46% and 0.48% for M1 and M2, resp.) and was significantly worse than model performance with actual tunings (p < 10 À30 and 10 À18 , resp.). Table 1 shows the beta coefficients and explained variances when factors were tested individually (regular font) or combined (bold face font) in the GLM. In both monkeys, ocular dominance (OD) was the least predictive of all factors tested and was nonsignificant for both correlation probability and correlation strength (b1, regular font). This factor dropped out upon initial screening. The weakness of OD was not simply because the arrays were aligned along ocular dominance columns. We confirmed this by recording from a third monkey in which the array was aligned across OD columns. In that monkey, we found no significant relationship between OD and correlation probability (p = 0.2) and a weak but significant relationship for correlation strength (R 2 = 0.008, p = 0.017, see Section 4).

Beta coefficients
Correlation probability was significantly associated with other receptive field properties including spatial frequency, phase, orientation, STRF, and full field color similarity, explaining 2.2-6.5% of the variance (b2-4, regular font). Correlation probability was more associated with STRF similarity, whereas correlation strength was more associated with spatial frequency/phase/orientation similarity. For STRF, STA-STA similarity alone predicted better than the combination of STA and STC, including models based on receptive field energy (1.5% and 2.3% for M1 and M2, resp.) and Procrustes transformation of STC components. Tuning similarity for spatial frequency or orientation alone was predictive of spike correlation, but phase alone was not (orientation: 7.4% and 4.1% for M1 and M2, resp.; spatial frequency: 3.6% and 2.8% for M1 and M2, resp.; phase was not significant: p = 0.26 M1, p = 0.09, M2). Color tuning was the best predictor of correlation strength (5.6-8.0% explained variance) and the second-best predictor of correlation probability (3.2-4.6%). Consistent with previous reports, spike correlation declined with distance (i.e. negative beta coefficients) when horizontal distance was tested as a single factor without cross-validation ( . ., N pairs, learned from the 'training' pairs only (no 'testing' pairs). x ik indicates tuning similarity of factor k (K = 5) for pair i. p i is defined as P(Z i = 1), where Z i is binary 0 or 1 depending on correlation significance. Because Z is difficult to visualize, we plotted a smoothed version of Z (dots, 'testing' pairs only) by convolving Z with a sliding window (width = 15). Although correlation probability differed across monkeys, the beta coefficients were similar across monkeys (Table 1). Goodness-of-fit measures were computed for each model by cross-validating across cells. (B and D) Model performance was assessed as area under ROC curve (thick curve: cross-validated across cells; thin solid curve: cross-validated across pairs). Axes indicate how the true and false positive rates for predicting correlation significance vary as a function of p, where ground truth is based on the correlation significance of independent data ('test' cells or pairs) not used to calculate p. In both cases, performance was significantly higher (p < 10 À10 , paired t-test, two-tailed, n = 20 cross-validations) than chance that is expected from columnar or laminar organization, modeled by vertical (within-penetration) and horizontal (within-row) shuffling (dashed and dotted curves, resp.). Generalized linear model beta coefficients based on tuning similarity and horizontal cortical distance (x 1 to x 5 ) were used to predict correlation probability (CP) and correlation strength (CS) measured during spontaneous activity. The beta coefficient of each factor is shown for the single-factor model (regular font) and for the combinedfactors model (boldface font, five factors). Model performance is cross-validated across cells and shown as percent explained variance (mean ± SEM) for the two monkeys. and K). With cross-validation, this relationship was weakly significant (p = 0.02 and 0.04 for M1 and M2) for correlation probability but was non-significant for correlation strength (p = 0.07 and 0.21 for M1 and M2, single-factor GLM, Table 1). The relationship became non-significant for multiple-factors GLM. We considered whether these results might be biased by uneven sampling with fewer pairs at larger distances. When the model was based on the same number of pairs at each distance (i.e. 64 pairs per distance for M1, 40 pairs per distance for M2), the beta coefficients for both distance and ocular dominance were weak and non-significant with cross-validation, for both correlation probability and correlation strength. Additionally, vertical distance and diagonal (contact-to-contact) distance were also nonsignificant factors in both monkeys (p > 0.1). Thus, the weakness of distance and ocular dominance in the model is not due to uneven sampling in distance.

Additivity/co-linearity of functional properties
Because some of these properties are related (e.g. STRF and SF/ phase/orientation), we wondered whether separating the factors individually or examining particular combinations of factors would lead to a different interpretation. We examined pair-wise correlation coefficients (R) between tuning similarities of these factors for each monkey ( Fig. 10A and B, for M1 and M2 respectively). The overall pattern showed that ocular dominance similarity was uncorrelated with other factors. STRF, color, and SF/phase/orientation similarity were positively correlated. STRF similarity and horizontal distance were negatively correlated ( Fig. 5E and F).
We also considered whether co-linearity may have resulted in non-stationary beta coefficients, e.g. between STRF and horizontal distance. However, varying these beta coefficients (varying b3 and b5 while keeping other beta coefficients constant) resulted in worse cross-validated performance (not shown). Thus, although the collinearity of STRF and horizontal distance resulted in some instability in the beta coefficients, it was still possible to obtain stable values for both b3 and b5.
Logically, collinearity should cause the total explained variance of combined factors ('Combined') to be less than the sum of the explained variances of individual factors ('Single'). Fig. 10C and D shows the summed explained variance for single factors, vs. the explained variance for the combined factors. As expected, collinearity reduced the additivity of these factors, resulting in lower (by 10-20%) explained variance for combined factors.
Next, we examined how different combinations of factors affected GLM performance (Fig. 11A-D). The additional contributions of individual factors varied depending on whether it was the only factor included in the model (slope of lines between 0 and 1 factor) or if it was the last factor added to the model (slope of lines between 4 and 5 factors). The ordering had little effect on the added contribution of most factors, but appeared to affect the contribution of horizontal distance (yellow). To quantify this additional contribution, we compared the performance of GLM when individual factors were left out of the combination (e.g. a 4-factor model), vs. when all 5 factors were combined. The additional contribution of each factor to performance (delta explained variance) is shown in Fig. 11E and 11F for M1 and M2, respectively, for both correlation probability (CP) and correlation strength (CS) (unpaired t-test, N = 20). Based on this test of additional contribution, horizontal cortical distance does not provide significant contributions beyond those provided by other factors, for both correlation probability and correlation strength (p > 0.1). We note that OD also does not provide an additional contribution, but it was already nonsignificant as a single factor upon cross-validation.
Finally, we considered informative predictors in our final model. During the model selection (see Section 2), only ocular dominance was removed during the initial screening due to its non-significance in the single factor model (p > 0.2). When we ran a combined factors model based on the remaining four factors, horizontal distance was non-significant (p > 0.05) and so was dropped. Reducing the combined factors model from 5 factors to the remaining 3 factors was also supported by likelihood ratio test between the two models ( Table 2). The reduced GLM contained only three remaining factors (x 2 -x 4 ), SF/phase/orientation, STRF, and color. Table 3 shows the beta coefficients, tuning similarity scores used as factor values (see tuning similarity histograms in Figs. 2G, 3G, and 4H), and total explained variance for this reduced model. Beta coefficients in the 3-factor model were also compared with those of the 5-factor model (unpaired t test). Notably, all beta coefficients were not significantly different across the two models, except for STRF (b3). This further indicates substantial collinearity between STRF and cortical distance. We suggest that these beta coefficient values may be useful for developing computational models of V1.

Discussion
Our results show that, for short distances up to 1.4 mm (about two V1 hypercolumns), tuning dissimilarity is better than cortical   Fig. 11. Contributions of combinations of GLM factors. Tree graphs show cell-wise cross-validated performance for all possible permutations of five factors (N = 5!) for correlation probability (A and B) and correlation strength (C and D) across the two monkeys. Starting from the left, the GLM's explained variance generally increases with increasing the number of factors, depending on the combination of factors included (e.g. blue line at arrowhead indicates how adding color improves the GLM performance beyond just including STRF). Note that in every case, adding STRF (red) improved performance, whereas adding horizontal distance resulted in negligible improvement after $3 factors are included. (E and F) The additional contribution of each factor was measured as the difference in explained variance when all 5 factors were combined in the model, vs. when that factor was dropped (4 factor model). ÃÃ p < 0.01, ÃÃÃ p < 0.001, n.s. p > 0.05. Although spike correlation does decline with cortical distance, this can be explained at short distances by the decline in tuning similarity, without requiring an additional factor for cortical distance itself. Thus, measuring combinations of tuning similarities may be necessary to model ''distance'' in the cortex, at least for short cortical distances. Also, previous studies showed that different tuning properties predict spike correlation, but it was unclear what was the relative contribution of these properties. Our model suggests that color (and luminance) tuning similarity is the strongest predictor, followed by spatiotemporal receptive field similarity. Surprisingly, ocular dominance was not a significant predictor beyond layer 4. We report these factor weights (Tables 1 and 3) and are placing this data in a public database (http://crcns.org/). This data is likely to be useful for constructing computational models of ensemble dynamics.

Why was horizontal distance non-significant?
Consistent with previous reports, raw spike correlations did decline with horizontal distance (Fig. 6H-I). The steep decline in raw correlations is associated with a decline in the incidence rate of fast correlations ( Fig. 6J and K), and a shallower decline in the strength of fast correlations ( Fig. 6L and M). This dependency on distance disappeared after a combination of tuning similarity factors was included. In summary, both the focus on fast correlations (by applying jitter correction) and the combination of multiple factors (Fig. 11) in a generalized linear model contributed to exclude short cortical distance as an explanatory factor.
Although this may appear to contradict a previous report that found that shuffle-subtracted spike correlations decreased with distance independent of orientation similarity (Das & Gilbert, 1999), we suggest that orientation similarity is only one of many tuning similarities (e.g. spatial frequency, phase, motion, ocular preference) that combine to determine spike correlation. In our factor permutation analysis (Fig. 11), distance is significant unless multiple factors are included.
An important facet of these findings is that they are based on arrays aligned along OD columns, and that anisotropy of horizontal projections may lead to different results for arrays aligned across OD columns. Unfortunately, currently available electrode arrays do not enable sufficiently dense and unbiased sampling of combinations of factors (ideally one would sample multiple depths per 80 lm horizontal distance, (Nauhaus et al., 2008)), and aligning our electrode arrays across OD columns tends to miss color blobs and oversample iso-orientation domains, resulting in fewer color-sensitive cells and stronger spike correlation strength at short distances due to the orthogonality of orientation and ocular dominance maps. These results are also limited to the short distances examined (up to 1.4 mm) and do not include factors that may apply at longer distances, e.g. long-distance patchy connections within V1. Longer distance interactions may be specific for particular properties (e.g. orientation, color) and may be associated with contextual effects (e.g. co-linearity, co-circularity) that are beyond the scope of this study.

Why was ocular dominance similarity non-significant?
Our ocular dominance (OD) results are also surprising, because some reports suggest that OD and spike correlations are linked. Three reasons may explain the apparent discrepancy of our OD results. First is that correlation probability (CP) and correlation strength (CS) are not the same as counting the number of significant pairs, which was often the basis of reports in physiological and anatomical studies (Hata et al., 1991;Yoshioka et al., 1996). OD is a significant factor when we simply analyze based on number of significant pairs, instead of CP or CS (see Fig. 8 and Results).
Second is that our analysis spans most cortical layers, which emphasizes common input that is possibly local, rather than thalamocortical inputs that arrive in layer 4. When we analyzed only layer 4 (both cells within layer 4), OD was clearly a factor for the number of significant pairs (Fig. 8B), CP (Fig. 8G), and CS (Fig. 8L).
Third is that we did not exclude binocular cells from the OD analysis. Previous studies excluded binocular cells (Hata et al., 1991;Ts'o, Gilbert, & Wiesel, 1986), which may have elevated the importance of OD or the contribution of layer 4.
We considered that the weakness of OD might be because the arrays were aligned along OD columns. In a third monkey, we aligned the array across alternating OD columns. We found no significant relationship between OD and correlation probability (p = 0.2) and a weak but significant relationship for correlation strength (r 2 = 0.008 cross-validated across cells, p = 0.017, singlefactor GLM). As with the first two monkeys, spike correlations depended more on OD when the analysis was restricted to layer 4 (Fig. S3).
Our results agree with a previous report (Chiu & Weliky, 2002), that correlated spontaneous activity in ferret V1 (P24-P29) does not solely reflect the pattern of segregated eye-specific LGN inputs. Like us, they found non-significant relationship between ocular dominance similarity and spontaneous correlations within the same contralateral band. In regions of alternating ocular dominance patches, they did find a weak link (r 2 = 0.06) between spontaneous correlations and ODSI. Ts'o, Gilbert, and Wiesel (1986) also reported weak but non-significant relationship for OD (p < 0.06). This weakness is consistent with previous anatomical reports ( Malach et al., 1993;Yoshioka et al., 1996) that horizontal axons target freely across ocular dominance columns and are not confined to columns from the same eye. Thus, local circuits appear to largely dilute the effect of eye-specific thalamocortical inputs. Total EV (%) 12.9 ± 0.6 8.5 ± 0.5 14.2 ± 1.0 13.2 ± 0.3 Reduced model, combining only 3 factors without OD and horizontal cortical distance. Model performance is cross-validated across cells and shown as percent explained variance (Total EV: mean ± SEM). Although the precise beta coefficients varied across monkeys, their relative magnitudes across tuning factors were consistent. Beta coefficients were mostly unchanged between 3 factor model vs. 5 factor model. Significant difference for STRF is likely due to collinearity with horizontal distance. * p < 0.05, unpaired t-test.
Additional factors may also explain the stronger link between OD and spike correlations. In previous studies, visual stimulation (i.e. during strong thalamocortical input) may have elevated the link between synchrony and ocular dominance similarity, especially if synchrony was measured during effectively monocular visual presentation (Kruger & Aiple, 1988), if the receptive fields from the two eyes were not precisely aligned (e.g. a separate bar shown to each eye, resulting in asynchronous input) (Ts'o, Gilbert, & Wiesel, 1986), or if OD was collinear with other factors like distance (Hata et al., 1991).

Color (and luminance) best predict correlation strength
Although STRF was the best predictor of correlation probability, color (including luminance) was the best predictor of correlation strength (Table 1 explained variance values). Why was color such a strong predictor, given that V1 is full of orientation selective cells? One possible explanation is that because color cells are less selective for orientation (Lu & Roe, 2008), the explanatory power of color is particularly strong for such pairs. This is consistent with several studies that reported a significant association between chromatic preference and spike correlations (Roe & Ts'o, 1999;Ts'o & Gilbert, 1988). It is possible that aligning our arrays along OD columns may also have increased the proportion of the population that was color-selective, due to the tendency of color blobs to align along OD columns (Lu & Roe, 2008). The predictive power of color similarity may be even stronger if color subfields were tested, instead of the full-field color patches we tested here.
Because our measure of full field color similarity is based on both hue and luminance, it is reasonable to ask whether most of it is due to hue or due to luminance? The answer is, ''both''. When assessed individually, neither hue nor luminance was a significant predictor, and it was only in combination that they were significant. We previously reported that spike correlation depended on luminance edge contrast similarity (Hung, Ramsden, & Roe, 2007), so there is good reason to suspect that luminance tuning may predict spike correlations. The results here suggest that hue and luminance act in combination to determine spike correlations, consistent with the integration of these cues in V1 (Clifford et al., 2003).

Generalizability of these results across cell tuning strength and animal state
Besides tuning similarity, spike correlation is also related to the degree to which the cells are tuned for particular factor (tuning strength). Like OD, the strength of color tuning is stronger in layer 4. However, we did not see an obvious difference in the number of significant pairs, CP, and CS for layer 4 vs. other layers (compare Fig. 8C,H, and M with Fig. 8D, I, and N). Instead, the main difference was a higher number of significant correlations between gray-gray neurons due to the prevalence of such neurons (Fig. 8C), but higher CP and CS between color-color pairs ( Fig. 8H and M). There also appears to be some dependency on cortical depth (Fig. 8E, J, and O), as well as on asymmetric properties (e.g. simple vs. complex, low vs. high spatial frequency tuning). We did not add these factors to the model because they were peripheral to the question of whether (short) cortical distance or tuning dissimilarity better explained spike correlations, and because they greatly increased the number of factors and added the complication of temporal sequence, while only modestly increasing performance. Overall, spike correlations are more complicated than we have described with this simple model, and other tuning properties are also important factors that are not included in our model.
Because spike correlations also depend on stimulus condition, we also examined spike correlation acquired from one monkey (M1) during the OD block. We chose the OD block because it activated most cells, including cells that preferred color, and also because we wanted to confirm the lack of relationship with OD during spontaneous activity. The results of the GLM were essentially the same, with only slightly lower explained variance in the OD block compared to spontaneous activity (logistic: 12.7%, logarithmic: 11.7%, cross-validated across cells). Specifically, OD remained non-significant, and horizontal distance was also non-significant when other factors were included. SF, phase, and orientation remained highly significant. Color tuning similarity remained significant for CS, but it was non-significant for CP, possibly due to different tuning for full-field color vs. the OD color gratings.
Finally, measurements of spike correlations may depend on the state of the animal. Spontaneous synchrony in barbiturate anesthesia, particularly at concentrations many times higher than we have used here (Contreras & Steriade, 1997), is qualitatively and quantitatively different from recordings in awake animals. Interactions under anesthesia reflect more the cumulative weight of common input and similarly tuned connections (Kenet et al., 2003;Tsodyks et al., 1999), whereas dynamic effects such as attention (Roelfsema et al., 1997;Singer, 1999) may emerge or dominate during natural vision. On the other hand, recordings in awake animals rely more on eye-segregated thalamocortical input and may be affected by the different microsaccades of the two eyes (Hermens & Walker, 2010), which may artificially elevate the dependence of synchrony on ocular dominance similarity. Spike synchrony and local field potentials are time-locked to saccades and microsaccades (Bosman et al., 2009;Ito et al., 2011;Rajkai et al., 2008). Recording under light anesthesia and muscle relaxation, as in this study, avoids this potential confound. In summary, combining multiple tuning properties measured during spontaneous activity can be useful to identify 'baseline' spike correlations to compare across studies, to avoid potential confounds, and to improve sensitivity for discovering what are the other non-similarity interactions (e.g. co-circularity) that are learned from natural image statistics.