## Abstract

Correlated spiking has been widely observed, but its impact on neural coding remains controversial. Correlation arising from comodulation of rates across neurons has been shown to vary with the firing rates of individual neurons. This translates into rate and correlation being equivalently tuned to the stimulus; under those conditions, correlated spiking does not provide information beyond that already available from individual neuron firing rates. Such correlations are irrelevant and can reduce coding efficiency by introducing redundancy. Using simulations and experiments in rat hippocampal neurons, we show here that pairs of neurons receiving correlated input also exhibit correlations arising from precise spike-time synchronization. Contrary to rate comodulation, spike-time synchronization is unaffected by firing rate, thus enabling synchrony- and rate-based coding to operate independently. The type of output correlation depends on whether intrinsic neuron properties promote integration or coincidence detection: “ideal” integrators (with spike generation sensitive to stimulus mean) exhibit rate comodulation, whereas ideal coincidence detectors (with spike generation sensitive to stimulus variance) exhibit precise spike-time synchronization. Pyramidal neurons are sensitive to both stimulus mean and variance, and thus exhibit both types of output correlation proportioned according to which operating mode is dominant. Our results explain how different types of correlations arise based on how individual neurons generate spikes, and why spike-time synchronization and rate comodulation can encode different stimulus properties. Our results also highlight the importance of neuronal properties for population-level coding insofar as neural networks can employ different coding schemes depending on the dominant operating mode of their constituent neurons.

## Introduction

Neurons in many brain areas exhibit correlated spiking but the role of those correlations remains controversial (Singer, 1993; Zohary et al., 1994; Engel et al., 1997; Gerstner et al., 1997; Shadlen and Movshon, 1999; Treisman, 1999; Salinas and Sejnowski, 2001; Palanca and DeAngelis, 2005; Averbeck et al., 2006; Schneidman et al., 2006; Wolfe et al., 2010). Noise correlations are generally thought to degrade coding efficiency (Averbeck et al., 2006) [with exceptions (Cafaro and Rieke, 2010)], but signal-dependent correlations could conceivably carry information. However, the feasibility of correlation-based coding has been called into question by the observation that output correlation varies with firing rate despite no change in input correlation (de la Rocha et al., 2007). If such a correlation–rate relationship always existed, input correlation could not be unambiguously decoded from output correlation, and transferred correlations would become meaningless (see Fig. 1). Importantly, correlations range from precise spike-time synchronization (on a millisecond timescale) to coarse rate comodulation (on a timescale up to seconds). We hypothesized that different types of correlation may differ fundamentally in how they are generated and what information they convey.

Propagation of correlated spiking depends on how individual neurons respond to correlated input (sensitivity to correlation) and whether groups of neurons respond with correlated output (transfer of correlation) such that postsynaptic neurons themselves receive correlated input (Abeles, 1991; Aertsen et al., 1996; Reyes, 2003). With respect to sensitivity to correlation, a critical factor is whether neurons operate as integrators or coincidence detectors: integrators respond to temporally dispersed inputs, whereas coincidence detectors respond selectively to rapid depolarization caused by temporally coincident (synchronous) inputs (Abeles, 1982; König et al., 1996). Operating mode (i.e., integration vs coincidence detection) reflects interplay between stimulus kinetics and spike threshold mechanism (see Results). With respect to the transfer of correlation, integrators spike repetitively at a rate proportional to their time-averaged input, whereas coincidence detectors respond to each suprathreshold input with an isolated spike—these spiking patterns are conducive to rate and temporal coding, respectively (Mainen and Sejnowski, 1995; König et al., 1996; Salinas and Sejnowski, 2001; Schreiber et al., 2004; Prescott et al., 2006; Prescott and Sejnowski, 2008; Tiesinga et al., 2008). Importantly, the temporal precision with which an individual neuron spikes should affect how synchronously a group of such neurons will spike to shared inputs, which bears directly on transfer of synchrony. Notably, neurons may be highly specialized for one operating mode, but most exhibit a context-dependent combination of modes (Maex et al., 2000; Rudolph and Destexhe, 2003; Prescott et al., 2006; Hong et al., 2008).

Using computer simulations and dynamic-clamp experiments, we identified different types of output correlation and investigated how transfer of each type of correlation depends on the intrinsic properties of cell pairs receiving correlated input. We show that pairs of “realistic” coincidence detectors exhibit rate comodulation and spike-time synchronization, whereas pairs of realistic integrators exhibit only rate comodulation. Synchrony, unlike rate comodulation, depends on spike-timing rather than rate. Rate- and synchrony-based coding are thus shown to operate independently.

## Materials and Methods

#### Stimulus preparation

We constructed stimulus waveforms using the same procedure as de la Rocha et al. (2007). The stimulus *I*(*t*) for a given trial was the linear summation of two Ornstein–Uhlenbeck (OU) processes (Uhlenbeck and Ornstein, 1930) described by the following:
where μ and σ are the mean and SD of the stimulus, and *c* is the input correlation (i.e., the fraction of fluctuating input shared between neurons) (see Fig. 1*A*). The common component ζ_{c}(*t*) was instantiated once and applied to all trials, whereas the independent component ζ_{i}(*t*) was randomly updated for each trial. Each OU process was formed by the following:
where ξ(*t*) is Gaussian white noise with zero mean and unit variance. Sampling rates were 10 and 5 kHz for experiments and simulations, respectively. *N*_{τ} = (2/τ)^{1/2} is a normalization constant that makes ζ(*t*) have unit variance. A correlation time τ = 5 ms was used unless otherwise indicated.

#### Model neurons and simulation procedures

Two conductance-based neuron models were used. We modeled the integrator as a Morris–Lecar (ML) model with type 1 excitability (Prescott et al., 2008a) and the coincidence detector as a Hodgkin–Huxley low-sodium (HHLS) model with type 3 excitability (Lundstrom et al., 2008). Equations for the ML model are as follows:
where
and *g*_{Na} = 20 mS/cm^{2}, *g*_{K} = 20 mS/cm^{2}, *g*_{L} = 2 mS/cm^{2}, φ = 0.15, *V*_{1} = −1.2 mV, *V*_{2} = 18 mV, *V*_{3} = 0 mV, and *V*_{4} = 10 mV. The membrane capacitance per area *C* was 2 μF/cm^{2} and the surface area was 100 μm^{2}.

Equations for the HHLS model are as follows:
with activation variables *m*, *n*, and *h* governed by the following:
where
*g*_{Na} = 41 mS/cm^{2}, *g*_{K} = 79 mS/cm^{2}, *g*_{L} = 0.3 mS/cm^{2}, and the membrane capacitance *C* = 1 μF/cm^{2} and the surface area was 100 μm^{2}.

The filter-and-threshold (FT) model consisted of three components: a linear filter to transform input to voltage, a voltage threshold, and an afterhyperpolarization (AHP). For the filter, we used the time derivative of a 15-ms-long Blackman filter, which was normalized to transform an input with variance 1 pA^{2} to an output with a variance 0.1 mV^{2}. The threshold was 1 mV and the AHP inserted for each spike had −0.5 mV amplitude and 30 ms decay time.

All simulations in conductance-based models were performed in NEURON (Hines and Carnevale, 1997). Simulations in the FT model were performed using custom Python scripts. All code will be made available on ModelDB. Each stimulus condition (*c*, μ, σ^{2}) was repeated 2–10 times for 30 min of simulated time. All (*c*, μ, σ^{2}) combinations used are summarized in Table 1. For Figure 5, the 100–200 simulation runs conducted for each model resulted in a very large amount of data, making the calculation of correlations (from up to ∼40,000 pairs) computationally challenging and the results difficult to present; therefore, we selected only pairs having the same stimulus mean, which resulted in 1000 and 2500 pairs for the integrator and coincidence detector, respectively. This criterion was not applied to experimental data (see Fig. 7) since there were fewer trials.

#### Slice preparation and electrophysiology

Experimental protocols were approved by the University of Pittsburgh Institutional Animal Care and Use Committee and have been previously described (Prescott et al., 2006). Briefly, adult male Sprague Dawley rats were anesthetized with intraperitoneal injection of sodium pentobarbital (50–75 mg/kg) and perfused intracardially with ice-cold oxygenated (95% O_{2} and 5% CO_{2}) sucrose-substituted artificial CSF (ACSF) containing the following (in mm): 252 sucrose, 2.5 KCl, 2 CaCl_{2}, 2 MgCl_{2}, 10 glucose, 26 NaHCO_{3}, 1.25 NaH_{2}PO_{4}, and 5 kynurenic acid. The brain was rapidly removed and sectioned coronally to give 300-μm-thick slices, which were kept in normal oxygenated ACSF (126 mm NaCl instead of sucrose and without kynurenic acid) at room temperature until recording.

Slices were transferred to a recording chamber constantly perfused with oxygenated (95% O_{2} and 5% CO_{2}) ACSF heated to 31 ± 1°C. Pyramidal neurons in the CA1 region of hippocampus were recorded in the whole-cell configuration with >70% series resistance compensation using an Axopatch 200B amplifier (Molecular Devices). Membrane potential (after correction for the liquid junction potential of 9 mV) was adjusted to −70 mV through tonic current injection. Intracellular recording solution contained the following (in mm): 125 KMeSO_{4}, 5 KCl, 10 HEPES, and 2 MgCl_{2}, 4 ATP (Sigma-Aldrich), 0.4 GTP (Sigma-Aldrich), as well as 0.1% Lucifer yellow; pH was adjusted to 7.2 with KOH. Pyramidal morphology was confirmed with epifluorescence after recording. All experiments were performed in 10 μm bicuculline methiodide (Sigma-Aldrich), 10 μm CNQX (6-cyano-7-nitroquinoxaline-2,3-dione) (Sigma-Aldrich), and 40 μm d-AP-5 (d-2-amino-5-phosphonovaleric acid) (Ascent Scientific) to block background synaptic activity.

Stimuli (see above) were injected into the recorded neurons through the patch pipette. To manipulate spike threshold mechanism (see Results), an artificial “shunt” conductance (*E*_{shunt} = −70 mV, *g*_{shunt} = 10 nS) was applied via dynamic clamp implemented with a Digidata 1200A ADC/DAC board (Molecular Devices) and DYNCLAMP2 software (Pinto et al., 2001) running on a dedicated processor as previously described (Prescott et al., 2006; Prescott and De Koninck, 2009); update rate was 10 kHz. Traces were low-pass filtered at 2 kHz and digitized at 10 kHz using a CED 1401 computer interface (Cambridge Electronic Design).

#### Reverse correlation analysis

For simulation and experimental data, we calculated spike-triggered averages (STAs) and covariance of stimuli (STC). The STA is simply the average of the set of stimuli that led to spikes subtracted from the mean of the prior stimulus distribution (i.e., the distribution of all stimuli independent of spiking output). Here, to remove the ambiguities caused by the temporal correlations, we used the fluctuating part of the unfiltered stimulus, *t*) instead of ζ(*t*) (see Eq. 2). Therefore, we have the following:
The time window for the STA was 200 ms before each spike, which captured most of the STA power. In a similar way, the STC and spike-triggered correlation of stimuli (STCor) *Q*(*t*,*t′*) are given by Bialek and de Ruyter Van Steveninck (1988) as follows:
where Cov_{spike} and Cov_{prior} are the covariance matrices of the spike-triggered and prior stimuli, respectively. The STA and STCor were used for predicting a cross-correlogram (CCG) in the first- and second-order by Equations 29 and 31, respectively (see below, Predicting CCGs from reverse correlation analysis). Predicted spike train covariance and correlation were computed from the CCG in the same way as the measured CCGs (see below).

#### Calculation of the measured CCG and correlation

We computed the CCGs of each spike train as follows. We started by building a spike train from the spike times with a Δ*t* = 1 ms time bin and computed the CCGs via the time-averaged unbiased empirical correlation function (Perkel et al., 1967): when a spike train for neuron *i* and *k*th repetition is *y _{i}*

_{,k}(

*t*), the cross-correlogram is given by the following: where

*N*

_{repeat}is the number of repetitions, and

*L*is the length of the spike trains and

*y*

_{i}_{,Nrepeat+1}=

*y*

_{i}_{,1}.

From Equation 10, we computed the correlation of the spike-counts with the time window of size *T* as follows:
ϴ* _{T}* is a window function giving ϴ

*(*

_{T}*t*) = 1 if 0 ≤

*t*≤

*T*and ϴ

*(*

_{T}*t*) = 0 otherwise. The shift correlator computes the covariance as follows: and the correlation coefficient is given by ρ

*=*

_{T}*C*/(

_{T}*C*

_{1}

*C*

_{2})

^{1/2}where the auto-covariance of each neuron

*C*(

_{i}*i*= 1,2) was calculated in the same way. Equations 10 and 12 are related as follows: where Φ

*(*

_{T}*t*) is a triangular window function such that Φ

*(*

_{T}*t*) = |

*T*−

*t*| if −

*T*≤

*t*≤

*T*and Φ

*(*

_{T}*t*) = 0 otherwise (Bair et al., 2001). Therefore,

*C*as well as ρ

_{T}*were computed from CCGs via Equation 13.*

_{T}In our analysis, we always computed the full CCGs and tried to analyze their behavior both at long and short timescales. The time window size *T* was not as important as in de la Rocha et al. (2007); therefore, in every case, we used *T* = 200 ms and dropped *T* from the notations.

#### Statistical analysis

We calculated the variance of the covariance, Var[*C*], by computing the bootstrap statistics of <*y*_{1}(0)*y*_{2}(τ)τ_{shuffle} in Equation 10; at each step, we constructed a resampled cross-correlogram <*y*_{1}(0)*y*_{2}(τ)τ_{resample} by random resampling from <*y*_{1}(0)*y*_{2}(τ)τ_{shuffle} and computed <*n*_{1}*n*_{2}τ_{resample} by Equation 12. We collected 400 <*n*_{1}*n*_{2}τ_{resample} with which the variance became stable. The Var[*C*] is taken as this resample variance and *C* is also recalibrated with the resample mean as follows:

For experimental data in Figure 7*B–D*, *t* scores were calculated from *C* and Var[*C*]. However, as for the simulation data, we could arbitrarily add more repetitions to suppress Var[*C*] down to sufficiently low level. For the analytic prediction, Var[*C*] could be computed in a similar way, but it was always insignificant (<1% of the covariance) for both experimental and simulation data.

In Figures 5 and 7, we evaluated the predictive power of the predicted *C* in terms of the coefficient of determination, *R*^{2}, defined by the following:
By definition, *R*^{2} = 1 signifies the perfect match, while the prediction has no correlation with the measured *C* when *R*^{2} = 0. In some cases, *R*^{2} < 0, and this signifies that the predictions in fact diverge away from the data with a different average.

#### Predicting CCGs from reverse correlation analysis

Here, we derive the first- and second-order prediction presented as Equations 29–32 in Results by using the reverse correlation analysis of the linear–nonlinear (LN) cascade model (Victor and Shapley, 1980; Meister and Berry, 1999). The LN model is composed of two stages: first, the stimulus *I*(*t*) was linear filtered by the relevant features of the model {ε_{α}} (α = 1, 2, …, *D*), and the probability to spike at the given time bin [*t*, *t* + Δ*t*], *P*(spike|**x**) where
Here, as above (see Reverse correlation analysis), we use the unfiltered stimulus *I*(*t*) are linearly related, STA and STCor have the same temporal correlation as the actual stimulus, and, consequently, the temporal correlation naturally shows up in the predicted CCGs; for example, compare predicted CCGs for τ = 5 and 50 ms in Figure 5*B*.

We follow the same derivation as in the study by de la Rocha et al. (2007): the common noise part δ*I* = *I*_{0} = *I*(*t*)|_{c}_{=0} and the comodulated part of the firing rate will be determined by Wiener series expansion in δ*I*.

##### First-order prediction (see Results, Eqs. 29, 30).

The firing rate change induced by a small perturbation *I*_{0} → *I*_{0} + δ*I* can be approximated as follows:
where δ*I*(*t*) has Gaussian statistics with a variance σ^{2} (Rieke et al., 1999; Hong et al., 2008). STA represents the spike-triggered average of stimuli. Then the correlation function of *q*(*t*) = *P*(spike at *t*) is as follows:
where <δ*I*_{1}(0)δ*I*_{2}(τ)τ = *c*σ_{1}σ_{2}δ(τ) if δ*I _{i}* =

We now define the firing rate *v*(*t*) = *q*(*t*)/Δ*t* and the predicted cross-correlogram for each pair as follows:
which is Equation 29 in Results. Furthermore, if we use the identity (Chialvo et al., 1997; Hong et al., 2008),
the firing rate covariance is given by the following:
where we changed the variables in the third line as *t″* = *t* + *t′*. The last line is equivalent to the original relationship between the correlation and gain (see Results, Eq. 30) (de la Rocha et al., 2007; Shea-Brown et al., 2008).

##### Second-order prediction (see Results, Eqs. 31, 32).

When we include the second-order term in Equation 17, we obtain the following (Rieke et al., 1999; Hong et al., 2008):
where *Q*(*t*_{1},*t*_{2}) is the STCor in Equation 9. The contribution to the CCG, Equation 31 in Results, is obtained in a straightforward way, as follows:
Note that we excluded the self-contraction to generate a proper Wiener series (Rieke et al., 1999).

In particular, when we have the same input condition, μ_{1} = μ_{2} = μ and σ_{1} = σ_{2} = σ, for the same type of neuron as in Figures 3 and 6, the peak height at *t* = 0 is given by the following:

Furthermore, we consider the case when the neuron is effectively well described by a single preferred spike-evoking stimulus feature (PSESF), ε(*t*), such as when the neuron functions almost as a pure integrator or a pure coincidence detector. In this case, *Q*(*t*_{1}*,t*_{2}) can be written as *Q*(*t*_{1},*t*_{2}) = *Y*(*t*_{1})*Y*(*t*_{2}), where *Y*(*t*) ∝ ε(*t*), and therefore
From ^{3}(∂*v*/∂σ)/*v* (Hong et al., 2008), the predicted peak height becomes
which is equivalent to Equation 32 in Results. Therefore, the firing rate gain with respect to stimulus variance contributes to the CCG peak height.

Note in the case of a single PSESF that the input mean and variance sensitivity of the firing rate are related as follows (Hong et al., 2008):
where
Therefore, when ε̄ ≠ 0, as in the pure integrator, the contribution of the peak to the correlation is suppressed compared with the first-order prediction, Equation 30. In the pure coincidence detector, on the other hand, ε̄ = 0, but firing rate ν is also independent of μ (see Fig. 6*B*), which can still make ∂ν/∂σ finite. Therefore, the predicted correlation in the pure coincidence detector [based on Eqs. 29, 30, and ν(μ) = constant] (Barreiro et al., 2010) will be profoundly modified by the quadratic-order approximation (see Results, Simulations in a phenomenological coincidence detector model).

## Results

When two or more neurons receive correlated (i.e., shared) input, their output spike trains should exhibit some correlation despite the effects of independent input (Fig. 1*A*). Spike train covariance *C* and the correlation coefficient ρ (i.e., *C* normalized by spike train variance) should, therefore, carry information about the input correlation *c*. However, de la Rocha et al. (2007) showed that the relationship between input and output correlation (denoted correlation susceptibility *S* = ρ/*c*) depends on the mean μ and variance σ^{2} of the input (Fig. 1*B*). This is important because stimulus-dependent changes in *S* prevent one-to-one mapping between input correlation and output correlation. Correlation-based coding is straightforward if *S* is independent of μ and σ^{2} (Fig. 1*C*, left), but it is compromised or requires a more complicated decoding scheme if *S* varies with μ and σ^{2} (Fig. 1*C*, right).

One way of understanding why this occurs is that, for a given input correlation, output correlation will vary depending on the sensitivity (i.e., gain) of the firing rate ν of each neuron with respect to the stimulus mean μ: if stimulus fluctuations occur within a steep region of the ν–μ curve, rate will fluctuate widely in each cell and the pair will exhibit large comodulated rate fluctuations, whereas stimulus fluctuations within a shallow region of the ν–μ curve will logically drive smaller comodulated rate fluctuations (Fig. 1*D*). Consequently, ρ “inherits” the same tuning as ν(μ), even when *c* is fixed; under these conditions, *c* cannot be unambiguously decoded from ρ without knowledge of μ. This line of reasoning triggered three concerns: (1) it applies to rate comodulation but not necessarily to spike-time synchronization, and thus it neglects one component of output correlation; (2) it implicitly assumes that input fluctuations are “noise” rather than “signal”; and (3) output rate and synchronization (i.e., correlation based on spike-time synchronization as opposed to rate comodulation) are liable to be tuned to different stimulus properties. This led us to our overall hypothesis that some cell types may encode signal-dependent fluctuations with precise spike-time synchronization and can do so independently of rate-based coding of other stimulus features.

The dependence of output rate ν on input parameters μ and σ differs fundamentally between cell types, as shown in Figure 2 for our conductance-based models: ν is principally sensitive to μ in the case of integrators, whereas it is also very sensitive to σ in the case of coincidence detectors (Higgs et al., 2006; Arsiero et al., 2007; Lundstrom et al., 2008). For simulations reported here, input was treated as a continuous stream rather than as discrete synaptic inputs (Destexhe et al., 2001); nonetheless, σ reflects coordinated fluctuations in presynaptic activity (Fellous et al., 2003), the temporal structure of which is reflected in the autocorrelation time τ (shorter τ implies more precise synchrony), while *c* specifies the proportion of inputs shared between two postsynaptic neurons. The first two parameters, σ and τ, affect the temporal precision of spiking in each postsynaptic neuron, while *c* affects correlation across neurons—all three parameters ultimately affect output synchrony.

The greater sensitivity of coincidence detectors to input synchrony relative to integrators (Fig. 2*A*,*B*) is a direct consequence of active membrane properties: activation of outward current (or inactivation of inward current) at perithreshold potentials helps ensure spike generation selectively in response to fast stimulus fluctuations (i.e., synchronous inputs), whereas perithreshold-activating inward current (or inactivating outward current) encourages repetitive spiking in response to constant or slow-changing input (Fourcaud-Trocmé et al., 2003; Svirskis et al., 2004; Higgs et al., 2006; Arsiero et al., 2007; Lundstrom et al., 2008; Prescott et al., 2008a). Differential sensitivity to input synchrony can be demonstrated most succinctly by contrasting which stimulus features preferentially elicit spikes in each cell type. We estimated the preferred spike-eliciting stimulus feature as the spike-triggered-averaged stimulus (STA) of the response of each cell to noisy input (see examples in Fig. 2*B*). The integrator exhibits a relatively broad, monophasic STA (Fig. 2*C*, left), whereas the coincidence detector exhibits a biphasic STA with a positive phase that is remarkably narrow (Fig. 2*C*, right).

The shape of the STA should, in theory, relate directly to the cross-correlation of output spiking given that the CCG is the overlap integral of the STA of each neuron, according to the first-order approximation in *c* (see Materials and Methods) (de la Rocha et al., 2007), as follows:
Hence, the broad monophasic integrator STA predicts a broad monophasic CCG for pairs of integrators (Fig. 2*D*, left), whereas the narrow biphasic coincidence detector STA predicts a narrow biphasic CCG for pairs of coincidence detectors (Fig. 2*D*, right). Using the relationship between the STA and firing rate, the covariance *C* of the output spike count is as follows (see Materials and Methods):
which is consistent with the study by de la Rocha et al. (2007) because sensitivity (or gain) of ν with respect to μ is the dominant factor determining the degree of correlated spiking given a certain degree of input correlation (Fig. 1*D*). Strong dependence of the relationship between ρ and *c* (i.e., correlation susceptibility *S*) on μ is not conducive to correlation-based coding (Fig. 1*C*, right). However, the STA does not always accurately represent the preferred spike-eliciting stimulus feature (this occurs when spike generation is sensitive to higher-order stimulus statistics), which invalidates predictions based on Equation 29 in certain cases (see below).

We therefore set out to identify (1) whether and how correlation susceptibility *S*(μ,σ) differs between integrators and coincidence detectors, (2) what the consequences of such differences are for correlation-based coding, and (3) precisely why the stimulus dependence of *S* differs between cell types.

### Simulations in conductance-based integrator and coincidence detector models

To compare the correlation susceptibility *S* (=ρ/*c*) of integrators and coincidence detectors, we conducted a series of numerical simulations using pairs of model neurons receiving a mix of correlated and independent input (Fig. 1*A*). For each neuron type, we varied mean μ under high- or low-variance σ^{2} conditions, and evaluated the output rate ν and correlation ρ (Fig. 3). As illustrated in Figure 2*A*, maximal ν differs between cell types and between σ conditions, which leads to unavoidable differences in the range of ν on each panel of Figure 3. Nevertheless, the ν(μ) curves (gray) are similarly shaped in all four conditions. If ρ simply inherited the same tuning as ν with respect to μ, then ρ(μ) curves should also be similarly shaped—clearly, they are not.

As expected for integrators, ρ/*c* increased to a maximum <1 as μ was increased in both the high and low σ cases (Fig. 3*A*, top). In coincidence detectors, however, a similar pattern was observed in the low σ case, but not in the high σ case, in which ρ/*c* decreased after reaching a peak (Fig. 3*A*, bottom). Figure 3*B* shows the same data replotted with ν on the *x*-axis. For coincidence detectors receiving high σ input, the correlation–rate relationship was strongly negative at high rates. The same coincidence detectors receiving low σ input exhibited little change in correlation across most of the (albeit narrow) range of firing rates, which implies that *S* is not significantly dependent on μ or ν in this cell type under these conditions. This property should be conducive to good correlation-based coding (see below).

To investigate how the stimulus dependence of *S* affects correlation-based coding by integrators and coincidence detectors, we measured ρ in response to different combinations of *c* and μ (Fig. 4*A*). As predicted, *S* was strongly dependent on μ in the case of integrators, which caused encoding of *c* by ρ to be ambiguous; in contrast, *S* was relatively unaffected by μ in the case of coincidence detectors, consistent with good correlation-based coding (compare Fig. 1*C*). Notably, Equation 29 failed to accurately predict ρ in the case of coincidence detectors (Fig. 4*A*, inset) despite having worked in the case of integrators (see below).

Small values of ρ observed among coincidence detectors stem from the shape of coincidence detector CCGs, which are narrow and biphasic, unlike integrator CCGs, which are broad and monophasic (Fig. 4*B*; compare prediction in Fig. 2*D*). Small values of ρ do not imply that synchrony transfer will fail, because downstream coincidence detectors prefer tall, narrow CCGs among their upstream neurons rather than short, broad CCGs, regardless of the absolute magnitude of ρ (S. Hong, S. A. Prescott, and E. De Schutter, unpublished observations). To summarize, correlation-based coding is viable despite small ρ as long as the relationship between ρ and *c* is insensitive to other stimulus features like μ. These data demonstrate the feasibility of correlation-based coding among coincidence detectors.

Next, we asked why *S*(μ,σ) differs between integrators and coincidence detectors. In all panels of Figure 3, ρ/*c* values predicted from Equation 29 were plotted for comparison with measured ρ/*c* values. The derivation of Equation 29 is based on linear response theory and requires several assumptions that are not strictly met in our simulations such as very small *c*. However, despite comparable stimulus parameters and identical methods used to calculate output correlation, Equation 29 predicted most of the output correlation observed for integrators but only a fraction of that observed for coincidence detectors. In the latter case, the degree of inaccuracy appeared to depend on σ. We reasoned that identifying why Equation 29 fails to accurately predict output correlation in coincidence detectors, especially under certain stimulus conditions, would help identify how correlated spiking in coincidence detectors differs from that in integrators.

To test the accuracy of the prediction, we plotted predicted covariance against measured covariance in Figure 5. If prediction by Equation 29 were consistently accurate (i.e., for all stimulus conditions), all data points would lie along the diagonal line. Figure 5*A* illustrates the consistent accuracy of the prediction for pairs of integrators versus its inconsistency for pairs of coincidence detectors. These data represent responses to a broad range of different μ and σ. The important observation is that the first-order prediction was reasonably accurate for all stimulus conditions in the case of integrators, whereas it was grossly inaccurate for many stimulus conditions in the case of coincidence detectors. Whereas Figure 3 illustrates stimulus conditions for which prediction by Equation 29 is good or bad, Figure 5 focuses on why Equation 29 sometimes fails to predict output correlation among coincidence detectors. By plotting covariance *C* rather than the correlation coefficient ρ, Figure 5*A* confirms that the prediction error seen for coincidence detectors is not attributable to autovariances in the spike train. Sample CCGs (Fig. 5*A*, right) reveal that the predicted CCG differs most from the measured CCG at the central peak; this is true for both integrators and coincidence detectors, but it translates into a larger discrepancy in predicting overall output correlation as the CCG gets narrower. Thus, these data qualitatively confirmed our starting predictions and identified a quantitative shortfall in the ability of Equation 29 to predict output correlation, especially in coincidence detectors.

In addition to reflecting the stimulus features to which each cell type is most sensitive, we reasoned that the width of the CCG also reflects the autocorrelation time of the correlated signal (which can be taken to reflect the temporal precision of spiking in presynaptic neurons) (see above). Therefore, we lengthened τ from 5 to 50 ms for simulations in Figure 5*B* based on the hypothesis that widening the CCG might reduce the prediction error. As expected, lengthening τ improved the prediction for both integrators and coincidence detectors but, whereas the CCG for the integrators was dramatically widened, the CCG for coincidence detectors remained narrow and the prediction still exhibited a sizeable error near the peak of the CCG. To determine whether prediction error near the peak of the CCG was sufficient to explain the discrepancy observed on the measured versus predicted *C* plots, we removed data within ±2 ms of the peaks of all coincidence detector CCGs and replotted the data for τ = 5 ms—this amounts to preferentially removing precisely synchronized spikes from the calculation of total output correlation. The result was a near-perfect prediction by Equation 29 (Fig. 5*C*), meaning the prediction error near the peak of narrow CCGs was sufficient to explain the failure of the prediction in coincidence detectors for all input conditions for which the original prediction was poor.

Next, we investigated (1) whether and how the original prediction might be improved so that the predicted CCG would more accurately match the measured CCG for coincidence detectors at the central peak, and (2) whether this would improve the measured versus predicted *C* plots (as expected given the results in Fig. 5*C*). We reasoned that if coincidence detectors are sensitive to stimulus variance (Fig. 2), then the prediction (Eq. 29) should take into account how stimulus variance affects spike generation. Therefore, we incorporated a second-order term into the prediction based on the following:
where *Q*(*t*,*t′*) is the spike-triggered correlation of the stimuli (STCor) (see Materials and Methods for derivation). Total output correlation should thus be predicted by the combination of Equations 29 and 31. For conditions in which input statistics and preferred spike-eliciting stimulus feature are equal, Equation 31 predicts the additional peak height as follows:
Using Equation 31, we find that the coincidence detector CCG was much better predicted at its central peak (Fig. 5*D*, right) and measured *C* was accurately predicted for a broad range of input conditions (Fig. 5*D*, left). Successfully modifying our quantitative prediction by incorporating Equation 31 argues that the difference in threshold mechanism is critical insofar as spike generation in integrators is preferentially sensitive to first-order stimulus statistics, whereas spike generation in coincidence detectors is also sensitive to second-order stimulus statistics (Fig. 2). On a more practical note, our ability to modify our quantitative prediction so that it works consistently (i.e., for all stimulus conditions) for integrators as well as for coincidence detectors confirms that we did not violate any fundamental assumption to the point of invalidating the prediction as a whole.

In principle, the second-order contribution (Eq. 31) can decrease much more quickly than the first-order prediction (Eq. 29) as the input correlation *c* becomes smaller. However, we found that the first-order prediction remained inconsistent for coincidence detectors even at relatively small values of *c* (Fig. 5*E*, left) and that including second-order terms significantly improved that prediction (Fig. 5*E*, right). In fact, measured *C* was closely matched by the second-order prediction across a broad range of *c*, and both were much larger than the first-order prediction except for very small *c* (Fig. 5*F*), proving that our result is not an artifact of large input correlations.

The nonlinear contribution to precise spike-time synchrony is not necessarily limited to second-order terms, and indeed the remaining systematic deviations from our second-order prediction in Figure 5, *D* and *E*, suggest that even higher-order terms contribute to synchrony, although that contribution is evidently quite small. Quantifying the *n*th order nonlinearity becomes more difficult for higher *n*, as the dimensionality of the *n*th moment of the spike-triggering stimuli rapidly increases, and was considered beyond the scope of the present study.

### Simulations in a phenomenological coincidence detector model

For the conductance-based models described above, the difference in threshold mechanism between integrators and coincidence detectors derives from distinct spike initiation dynamics (Prescott et al., 2008a). Based on results from phenomenological models without “dynamics” (de la Rocha et al., 2007), one might suspect that the specific dynamics of the cell model are not vital for the correlation–rate relationship. The type of threshold and the dynamical mechanism responsible for “thresholding” are inextricably linked, but, nonetheless, to test whether the type of threshold is critical rather than the dynamics per se, we constructed a phenomenological coincidence detector model comprising a filter and threshold in which the filter is biphasic, like the STA in our conductance-based coincidence detector model (Fig. 6*A*; compare Fig. 2*C*). In this dynamics-free model, ν was completely independent of μ but varied with σ (Fig. 6*B*; compare Fig. 2*A*, right)—in effect, this model constitutes an “ideal” coincidence detector. Despite the first-order prediction (i.e., Eq. 29) that there should be no output correlation (given that ∂ν/∂μ = 0) (Barreiro et al., 2010), pairs of ideal coincidence detectors receiving correlated input did indeed exhibit correlated spiking (Fig. 6*C*, left). Moreover, sample CCGs (Fig. 6*C*, right) were comparable in shape to those of the conductance-based coincidence detector model (compare Figs. 4, 5). Furthermore, the same prediction error was observed near the peak of the narrow CCGs when applying Equation 29, but our quantitative prediction was near-perfect when second-order terms (Eq. 31) were included. Overall, these results argue that the type of threshold used by a pair of neurons receiving correlated input will impact the correlation in their output.

To summarize results up to this point, ideal coincidence detectors (exemplified by our FT model) exhibit correlations arising solely from spike-time synchronization. In contrast, more realistic coincidence detectors (exemplified by our conductance-based HHLS model) exhibit a mixture of spike-time synchronization and rate comodulation, whereas realistic integrators (exemplified by our conductance-based ML model) exhibit mostly rate comodulation. The timescales of these two types of output correlation differ but nonetheless overlap. More importantly, the two types of output correlation exhibit fundamentally different sensitivities to firing rate: rate comodulation varies with firing rate, whereas spike-time synchronization does not. We predicted that real neurons should exhibit rate comodulation and spike-time synchronization if those neurons operate at least partially as coincidence detectors and, moreover, that the predominance of each type of output correlation will depend on the balance of operating modes.

### Experiments in CA1 pyramidal neurons made to behave preferentially as integrators or coincidence detectors via manipulation of threshold mechanism by dynamic clamp

Integrators and coincidence detectors are found throughout the nervous system but tend to exhibit specializations beyond threshold mechanism. Because those additional specializations could confound our comparison of correlation susceptibility, we chose to compare correlation susceptibility in a single type of neuron whose threshold mechanism was experimentally manipulated such that the neuron behaves preferentially as an integrator or coincidence detector.

Pyramidal neurons, including those in the CA1 region of hippocampus, display the hallmarks of integrators when recorded in acute brain slices but behave more like coincidence detectors upon introduction of a virtual leak conductance by dynamic clamp (Fig. 7*A*), consistent with a predicted switch in threshold mechanism (Prescott et al., 2006, 2008b). Thus, by manipulating the membrane properties of CA1 pyramidal neurons, we were able to compare correlation susceptibility in a single cell type operating preferentially in one or the other mode—the shift in operating mode is quantitative, not absolute. In the interests of comparing experimental and simulation data [and to compare with past studies (de la Rocha et al. (2007)], we applied the same stimulation paradigm used in simulations (i.e., noisy current injection) together with the simplest dynamic-clamp manipulation capable of switching the threshold mechanism (i.e., constant leak conductance).

Compared with STAs in the integrator-mode (Fig. 7*A*, left), STAs in the coincidence detector-mode (Fig. 7*A*, right) had a higher peak height, and were therefore steeper than integrator STAs, but the negative phase was far less prominent than in STAs from our coincidence detectors models (compare Fig. 2*C*). This is likely due to spike initiation in real neurons being influenced by more currents with gating variables spanning a broader range of timescales than were included in our minimal conductance-based models. In any case, as we would predict, the steep weakly biphasic STAs are consistent with the narrow weakly biphasic CCGs observed for pyramidal neurons in the coincidence detector-mode (see below). Most importantly for our purposes, firing rate gain with respect to μ, which is proportional to how biphasic the STA is (Eq. 20), was significantly lower in the coincidence detector-mode compared with the default integrator-mode (0.14 ± 0.041 vs 0.05 ± 0.016 Hz/pA; *p* < 0.01, *t* test), whereas firing rate gain with respect to σ was similar between modes (0.034 ± 0.022 vs 0.022 ± 0.011 Hz/pA), which means ∂ν/∂σ was, relative to ∂ν/∂μ, higher in the coincidence detector-mode (0.25 ± 0.16 vs 0.50 ± 0.29; *p* < 0.01, *t* test). This is consistent with the differential response properties reported for the conductance-based models in Figure 2*A*. Like for simulation data in Figure 5, we plotted our experimental data as measured versus predicted *C* based on the first-order prediction described in Equation 29 (Fig. 7*B*, top). As expected, the prediction error was much greater for the coincidence detector-mode than for the integrator-mode, which was also evident in the sample CCGs (Fig. 7*B*, bottom). Consistent with simulations (compare Figs. 4, 5, CCGs), the CCG for the coincidence detector-mode was much narrower than for the integrator-mode, and the first-order prediction deviated from the measured CCG primarily at its peak (Fig. 7*B*, bottom). As for simulation data, we hypothesized that removing a narrow region around the peak of the CCG (where the actual and predicted CCGs differ most) would dramatically improve the prediction as visualized on the measured versus predicted *C* plots, which indeed it did (Fig. 7*C*). Similarly, including the second-order terms described by Equation 31 improved our prediction, especially for the coincidence detector-mode (Fig. 7*D*). These experiments demonstrate that pyramidal neurons receiving shared input can exhibit output correlation in excess of that predicted from firing rate comodulation, and that this spike-time synchronization is greater when the neurons behave more like coincidence detectors.

## Discussion

Through simulations and experiments, we have shown that pairs of neurons receiving correlated input can exhibit spiking that is correlated in different ways, and that the overall contribution of each type of output correlation depends on intrinsic cellular properties. Ideal integrators exhibit output correlations comprised entirely of rate comodulation. This makes sense given that individual integrators use rate encoding in which spiking is dictated by the mean (i.e., a first-order statistic) of the input; therefore, input can be reconstructed entirely by applying a first-order filter (i.e., the STA) to the spike train. Conversely, ideal coincidence detectors exhibit output correlations comprised entirely of spike-time synchronization (Fig. 6). However, realistic coincidence detectors exhibit output correlations comprised partly of spike-time synchronization and partly of rate comodulation (Fig. 5). This too makes sense given that individual coincidence detectors use temporal encoding in which spiking is sensitive to the variance (i.e., a second-order statistic) of the input; therefore, reconstructing the input from the spike train requires inclusion of first- and second-order filters (Theunissen and Miller, 1995). Real neurons and realistic conductance-based models are never pure integrators or coincidence detectors. We have juxtaposed the two operating modes for didactic purposes, but, ultimately, operating mode is a continuum representing the interplay between stimulus kinetics and neural dynamics. Neural dynamics—most notably spike generation—differ between neurons, can be modulated within a given neuron, and directly impact how a neuronal ensemble encodes and transmits information.

Some cell types are optimized for integration or coincidence detection; for example, many neurons early in the auditory pathway (e.g., cochlear nucleus and superior olive) are exquisite coincidence detectors (Manis and Marx, 1991; Cao et al., 2007; Mathews et al., 2010), whereas certain neurons in the entorhinal cortex are near-perfect integrators (Egorov et al., 2002). In addition to morphological specializations (Agmon-Snir et al., 1998), a broad range of ion currents can influence operating mode, the common feature being that such currents activate or inactivate at voltages near threshold, or that they impact the gating of other perithreshold currents, for example, by shifting voltage threshold (Prescott et al., 2006, 2008a). Whether pyramidal neurons in the neocortex and hippocampus function as integrators or coincidence detectors has been controversial (Abeles, 1982; Softky and Koch, 1993; Shadlen and Newsome, 1998). The answer depends not only on intrinsic cell properties but also on stimulus conditions (Rudolph and Destexhe, 2003) and other external factors like background synaptic activity that influence membrane conductance (Destexhe et al., 2003). Shunting and adaptation, especially in combination, can encourage coincidence detection in pyramidal neurons that might otherwise behave preferentially as integrators (Prescott et al., 2006, 2008b). The adaptation current *I*_{M} has been shown to encourage coincidence detection in hippocampal and neocortical pyramidal neurons (Hu et al., 2007; Guan et al., 2011) and is itself subject to endogenous modulators and to several drugs (for review, see Brown and Passmore, 2009). This supports the view that pyramidal neurons (and presumably other cell types through adaptation by *I*_{M} and via other mechanisms) can shift between operating modes. By not operating at one or the other extreme, pyramidal neurons can flexibly use rate- and/or synchrony-based coding (see below).

The coexistence of rate- and synchrony-based coding is consistent with recent modeling work showing that differently encoded information can be simultaneously propagated through feedforward networks depending on network properties (Kremkow et al., 2010; Kumar et al., 2010). Our results emphasize the importance of cellular properties for exactly the same issues. Network and cellular properties can be intimately related (see above regarding synaptic input and membrane conductance) and may interact nonlinearly such that forms of modulation, which individually have weak effects, combine to produce powerful gating mechanisms that switch a network between propagating synchronous or asynchronous activity (see below). But it is not necessarily the case that synchronous and asynchronous spikes act independently; for example, in spiny stellate cells receiving thalamocortical input via a limited number of weak synapses, asynchronous background cortical input may be crucial for setting the membrane potential of stellate cells close enough to spike threshold that synchronous thalamic inputs can elicit spikes (Douglas and Martin, 2007; da Costa and Martin, 2011). In that example, asynchronous input may facilitate the propagation of synchronous spiking and, furthermore, might act as a continuously variable gain-control mechanism rather than as a simple on–off switch.

An important related issue is correlation between excitatory and inhibitory input. This has been proposed as a mechanism to regulate information transmission (Salinas and Sejnowski, 2001; Kremkow et al., 2010), amplify external inputs (Murphy and Miller, 2009), and enhance response fidelity (Wehr and Zador, 2003; Cafaro and Rieke, 2010). Such correlations may exist with delays of only a few milliseconds, likely because of feedforward inhibition ensuring a narrow integration window (Wehr and Zador, 2003; Higley and Contreras, 2006). Under those conditions, presynaptic spikes synchronized with millisecond precision are required to reliably evoke responses in postsynaptic neurons.

Information encoded by rate or synchrony must be reliably transmitted (when allowed by gating mechanisms) and eventually decoded. Unambiguous decoding requires independent rate- and correlation-based coding (Fig. 1). Recent findings (de la Rocha et al., 2007) have cast doubt on this independence. We concur that a correlation–rate relationship compromises correlation-based coding based on rate comodulation, but we demonstrate here that output correlation and rate can be correlated without interfering with synchrony-based coding. Counting spikes neglects the information carried by spike timing. The absolute timing of spikes in a single presynaptic cell has little impact on postsynaptic activation, but the relative timing of spikes across multiple presynaptic cells is crucial if the postsynaptic cell operates as a coincidence detector; by extension, the absolute timing of synchronized volleys likely carries information about the original stimulus.

It follows that integrators exhibit the correlation–rate relationship described by de la Rocha et al. (2007), but that same relationship, although it can be observed in coincidence detectors (Barreiro et al., 2010), applies to only a fraction (i.e., to only one component) of their output correlation. Notably, de la Rocha et al. used computer models in which spike generation depended on mean input (i.e., integration) and the 1-s-long stimuli used in their experiments were sufficiently short that slow adaptation was minimal, thus favoring spike generation based on mean input (Prescott and Sejnowski, 2008). In addition to testing “integrator” models, we tested conductance-based and phenomenological models in which spike generation depends on input variance (i.e., coincidence detector models). Our experiments were conducted with 300-s-long stimuli to favor adaptation, our rationale being that neurons *in vivo* exist in a chronically depolarized and shunted state because of synaptic bombardment (Destexhe et al., 2003). By testing a single cell type under different “virtual network conditions” rather than testing different cell types representing each extreme of the continuum between integration and coincidence detection, our results show that the neuronal operating mode can be shifted along that continuum and that this shift in operating mode adjusts the balance of rate- and synchrony-based coding.

These results demonstrate the importance of cell-level properties for network-level coding. There has been a longstanding bias in the network modeling community to focus on synaptic weights and network architecture, with far less emphasis put on cellular properties. Here, we have shown what happens when neurons do not “integrate and fire.” Similarly, Burak et al. (2009) and Barreiro et al. (2010) recently showed differences in output correlation depending on the type of model used. Effects of synaptic kinetics, background noise, and input spike train statistics on output correlation have also been documented (Maex et al., 2000; Tetzlaff et al., 2008; Ostojic et al., 2009). Tchumatchenko et al. (2010) have shown specifically that output correlations can be rate dependent or independent according to input conditions. Our results emphasize how both types of correlations can coexist and that this depends on input conditions and (how that input is encoded given) single cell firing properties. Indeed, whereas simplified models tend to favor pure integration or pure coincidence detection, real neurons and realistic conductance-based models exhibit a context-dependent mixture of operating modes. Neglecting the richness of cell-level coding properties will surely translate into underestimation of network-level coding possibilities.

In conclusion, different types of neurons or even the same neuron operating under different conditions can be differentially sensitive to first- and second-order stimulus statistics, namely mean and variance. If we do not assume that fluctuations are simply noise, and if a neuron is sensitive to those fluctuations (which is true of coincidence detectors), then output spiking should carry information about the signal variance. Two such neurons receiving correlated input may carry information about signal variance in their precisely synchronized spiking. Even if correlated with firing rate (indeed, such output correlations might arise from correlations between features of the original input), spike-time synchronization is separate from rate comodulation, which, by definition, is directly linked to the firing rate that is tuned to the signal mean. Both forms of correlation may coexist and operate independently. Thus, the richness of single-neuron coding abilities translates into even richer multineuron coding possibilities.

## Footnotes

This work was support by the Okinawa Institute of Science and Technology Promotion Corporation (S.H., E.D.S.) and NIH Grant NS076706 (S.A.P.). S.A.P. is also a Rita Allen Foundation Scholar in Pain and the 53rd Mallinckrodt Scholar.

The authors declare no competing financial interests.

- Correspondence should be addressed to Sungho Hong, Computational Neuroscience Unit, Okinawa Institute of Science and Technology, 7542 Onna, Onna-son, Okinawa 904-0411, Japan. shhong{at}oist.jp