# Intrinsic heterogeneity in oscillatory dynamics limits correlation-induced neural synchronization

## Abstract

Synchronous neural oscillations are found throughout the brain and are thought to contribute to neural coding and the propagation of activity. Several proposed mechanisms of synchronization have gained support through combined theoretical and experimental investigation, including mechanisms based on coupling and correlated input. Here, we ask how correlation-induced synchrony is affected by physiological heterogeneity across neurons. To address this question, we examined cell-to-cell differences in phase-response curves (PRCs), which characterize the response of periodically firing neurons to weak perturbations. Using acute slice electrophysiology, we measured PRCs across a single class of principal neurons capable of sensory-evoked oscillations in vivo: the olfactory bulb mitral cells (MCs). Periodically firing MCs displayed a broad range of PRCs, each of which was well fit by a simple three-parameter model. MCs also displayed differences in firing rate-current relationships and in preferred firing rate ranges. Both the observed PRC heterogeneity and moderate firing rate differences (∼10 Hz) separately reduced the maximum correlation-induced synchrony between MCs by up to 25–30%. Simulations further demonstrated that these components of heterogeneity alone were sufficient to account for the difference in synchronization among heterogeneous vs. homogeneous populations in vitro. Within this simulation framework, independent modulation of specific PRC features additionally revealed which aspects of PRC heterogeneity most strongly impact correlation-induced synchronization. Finally, we demonstrated good agreement of novel mathematical theory with our experimental and simulation results, providing a theoretical basis for the influence of heterogeneity on correlation-induced neural synchronization.

- phase-response curve
- stochastic synchrony
- neural oscillators
- biophysical diversity
- olfaction

synchronized oscillatory firing is a prominent feature of neuronal activity (Buzsáki and Draguhn 2004). Oscillatory synchrony can be task and state dependent (Salinas and Sejnowski 2001), is believed to critically contribute to neural coding (Panzeri et al. 2010), and is altered in several disorders (Uhlhaas and Singer 2006). Considerable progress has been made in elucidating mechanisms of oscillatory synchronization, particularly through combined theoretical and experimental approaches (Wang 2010). Direct coupling of inhibitory neurons or mixed populations of excitatory and inhibitory neurons by electrical and chemical synapses can synchronize oscillating neurons (Traub et al. 2004), and physiological heterogeneity can destabilize such coupling-induced synchronization (Kopell and Ermentrout 2002). Alternatively, fluctuating correlated input, such as that from common synaptic partners, can synchronize coupled or uncoupled oscillating neurons (Galán et al. 2006). How physiological heterogeneity in biophysical properties and firing rates affects correlation-induced synchronization is currently unknown.

The stereotyped circuitry and sensory-evoked oscillations of the rodent olfactory bulb (OB) affords an excellent model for studying the impact of heterogeneity on oscillatory synchronization (Kay et al. 2009; Padmanabhan and Urban 2010). Sensory-evoked excitation and asynchronous lateral inhibition drive mitral cells (MCs), the OB principal neurons, with correlated fluctuating input. Oscillating MCs are thus subject to both inhibitory coupling and noisy correlated input (Friedman and Strowbridge 2003; Galán et al. 2006; Lagier et al. 2004; Schoppa 2006). Moreover, intrinsic diversity within MCs, as measured through spike-triggered average (STA) heterogeneity, can significantly impact MC spike train correlations and information content (Padmanabhan and Urban 2010).

The degree to which oscillating neurons synchronize depends on their response to transient inputs, as captured by the phase-response curve (PRC) (Smeal et al. 2010). The PRC describes the response (phase advance or delay) of an oscillator to small transient inputs as a function of the phase at which the input arrives. Neurons in which the phase of spiking is always advanced by positive inputs have type I PRCs. Neurons in which positive inputs can either advance or delay spike phase have type II PRCs (Hansel et al. 1995). Previously, we have shown that homogeneous neurons with identical type II PRCs synchronize more strongly when driven by correlated inputs than do neurons with type I PRCs (Abouzeid and Ermentrout 2009; Galán et al. 2007; Marella and Ermentrout 2008). We hypothesize that, in turn, physiological levels of PRC heterogeneity significantly limit correlation-induced synchronization.

In this article, we present new data and analysis quantifying the influence of physiological heterogeneity on correlation-induced synchrony. In the first half of our study, we quantify the intrinsic biophysical diversity within oscillating MCs in vitro using acute slice electrophysiology. In the second half, we analyze how this observed heterogeneity, as well as systematic PRC heterogeneity in general, influences correlation-induced synchrony. We find that MCs firing in a roughly periodic manner (similar to odor-evoked oscillations observed in vivo) exhibit a broad range of type II PRCs that are surprisingly well characterized by a generalized three-parameter model. MCs also vary in their firing rate-current (FI) relationships and evoked periodic firing rates. Intrinsic biophysical diversity within this single class of principal neurons thus manifests as substantial and largely independent PRC and firing rate heterogeneity. Using within-cell comparisons to define the maximum synchrony possible among homogeneous MCs, we find that PRC differences between heterogeneous MCs reduce output synchrony by up to 30%. This is equivalent to the reduction in output correlation caused by reducing input correlation by 0.1. In turn, moderate firing rate differences (∼10 Hz) impose up to a comparable 25% reduction in maximum output synchrony. Combined, PRC and firing rate heterogeneity reduce maximum output synchrony by up to 40%, equivalent to reducing input correlation by as much as 0.2. Simulations further demonstrate that these two factors (physiological levels of PRC and firing rate heterogeneity) completely account for the difference in correlation-induced synchrony among heterogeneous vs. homogeneous populations in vitro. Independent modulation of PRC features reveals that correlation-induced synchrony among neurons with small-amplitude type II PRCs is more susceptible to PRC heterogeneity than synchrony among neurons with large-amplitude type II PRCs. Finally, we extend our previous analytical theory describing correlation-induced synchronization of uncoupled oscillators (Marella and Ermentrout 2008) to account for intrinsic oscillator heterogeneities. Predictions made with this generalized formulation closely match our experimental results, providing a theoretical context for how intrinsic diversity impacts correlation-induced synchronization of uncoupled oscillators.

## MATERIALS AND METHODS

### Physiological Experiments

#### Slice preparation.

Coronal OB slices (310–320 μm thick) were prepared from postnatal *day 12–18* mice of both sexes, as described previously (Giridhar et al. 2011). Briefly, C57BL/6 mice were anesthetized with isoflurane and decapitated into ice-cold oxygenated Ringer's solution containing (in mM) 125 NaCl, 25 glucose, 2.5 KCl, 25 NaHCO_{3}, 1.25 NaH_{2}PO_{4}, 1 MgCl_{2}, and 2.5 CaCl_{2}. OBs were isolated, and slices were obtained using a vibratome (VT1200S; Leica). Slices were removed to 37°C oxygenated Ringer solution for ∼15 min before recovering at room temperature for at least 30 min before electrophysiological recording. All experiments were completed in compliance with the guidelines established by the Institutional Animal Care and Use Committee of Carnegie Mellon University.

#### Electrophysiology.

Slices were continuously superfused with oxygenated Ringer's solution warmed to 37°C. MCs were identified by morphology and laminar position using infrared differential interference contrast (IR-DIC) microscopy. Unpaired whole cell current-clamp recordings were made from individual MCs using electrodes filled with (in mM) 120 K-gluconate, 2 KCl, 10 HEPES, 10 Na-phosphocreatine, 4 Mg-ATP, 0.3 Na_{3}GTP, and 1 EGTA. Synaptic transmission was not blocked, based on our previous observation that intrinsic MC biophysical diversity, as measured by *1*) STA heterogeneity, *2*) firing response across MCs to an identical input current waveform, and *3*) heterogeneity in FI relationships, is independent of intact synaptic transmission in unpaired recordings (Padmanabhan and Urban 2010). The insensitivity of STAs to synaptic blockade in particular is critical to our current results, because we use the STA to calculate the PRC for each oscillating MC (see below). Others have also noted that ongoing synaptic transmission in acute slice preparations does not strongly impact the accuracy of PRC estimation (Netoff et al. 2012). Data were low-pass filtered at 4 kHz and digitized at 10 kHz using multiclamp 700A and 700B amplifiers (Molecular Devices) and an ITC-18 acquisition board (Instrutech) controlled by custom software written in IGOR Pro (Wavemetrics).

#### Stimulus generation.

To study synchronization of oscillating neurons, we drove MCs to fire in a roughly periodic manner by injecting constant step currents (150–700 pA). These step currents were overlaid with frozen colored noise (σ = 15–20 pA) to *1*) more realistically mimic physiological synaptic input; *2*) estimate the STA, a biophysical metric unique to each neuron (Padmanabhan and Urban 2010) and directly related to the PRC for interspike interval (ISI) coefficient of variations (CV_{ISI}) <0.4 (Ermentrout et al. 2007); and *3*) examine the ability of small-amplitude, correlated fluctuating input to synchronize spiking (i.e., correlation-induced, “noise-induced,” or “stochastic” synchrony). To study correlation-induced synchrony, pairs of white noise inputs, ξ_{1,2}(*t*), with varying levels of correlation were produced by mixing common and independent white noise sources via ξ_{j}(*t*) = *C*_{in}ξ_{C}(*t*) + (1 − *C*_{in})η_{j}(*t*), where ξ_{C}(*t*) is a source of common noise and η_{1,2} are two sources of independent noise. *C*_{in} governs the correlation between the resulting pair of noise inputs, with *C*_{in} = 0 corresponding to no input correlation and *C*_{in} = 1 corresponding to perfect input correlation. ξ_{1,2}(*t*) were then normalized by 1 − 2*C*_{in} + 2*C*_{in}^{2} to maintain equal spectral properties and convolved with an alpha function: α(*t*) = (*t*/τ)exp(−*t*/τ) with τ = 3 ms. Constant step current injections were then overlaid with these correlated colored noise inputs.

### Data Analysis

#### PRC estimation.

To estimate the PRC for each in vitro MC, we took advantage of the direct relationship between the STA and PRC of a neuron firing in a roughly periodic manner (CV_{ISI} < 0.4) (Ermentrout et al. 2007; Torben-Nielsen et al. 2010). STAs were determined first, as follows. Each MC was injected with a series of noisy step current injections (see above), each 6 s in duration with 10 s separating each stimulus presentation. Binary strings of spike times were extracted from each resulting voltage trace. The initial 1 s of each trace was discarded to minimize effects of spike-frequency adaptation, and the remaining spikes were used to calculate CV_{ISI}. Under these conditions and in response to the described inputs, ∼80% of MCs exhibited noisy periodic firing (CV_{ISI} < 0.4) and were used for PRC and synchrony analyses. The STA of each MC firing in this noisy periodic regime was calculated as the mean injected current (with step current subtracted) preceding each spike by the natural period of the cell, defined as the mean ISI. STAs were calculated over 805–4,090 spikes (μ = 1,951 spikes; σ = 690 spikes), which well exceeds the few hundred spikes necessary for accurate PRC estimation using the STA method (Torben-Nielsen et al. 2010). To confirm that a sufficient number of spikes were used in each STA extraction, we compared the final STA with “trial-subsampled” STAs extracted from randomly sampling half of the total number of spikes (Padmanabhan and Urban 2010). To confirm consistency of the extracted STA across the duration of the recording, we compared the final STA with “time-subsampled” STAs extracted from the first and second halves of all spikes (Padmanabhan and Urban 2010). All final STAs closely agreed with their respective trial- and time-subsampled STAs (*R*^{2} > 0.95 for all cells), confirming accurate and consistent STA extraction.

The PRC of each oscillating MC was then derived from the STA using
*T* is the natural period (mean ISI), σ is the standard deviation of the injected current fluctuations, and *R*_{in} is the input resistance (Ermentrout et al. 2007). Using this formula, time *t* is mapped into phase θ in terms of radians by the conversion factor 2π/*T* to reflect the periodic nature of activity. This method yields a deterministic PRC for each in vitro MC; phase-dependent PRC variability was not examined but should not significantly influence the accuracy of our PRC estimate (Ermentrout et al. 2011). To systematically evaluate differences in PRCs across neurons, each PRC was fit to the generalized phenomenological model
*B* parameter governing the node or zero point of the curve. An amplitude component with parameter *A* describes the reactivity of the cell to phasic inputs, with higher values of *A* yielding greater advances and delays in spike timing for a given input. Finally, an exponential component governs the balance between phase advance and delay regimes, with high values of parameter *C* corresponding to little or no phase delays (i.e., closer to type I PRCs). This model provided highly accurate fits to all PRCs recorded in vitro (*R*^{2} = 0.9979 ± 0.0013, mean ± SD), facilitating quantitative analysis of PRC heterogeneity. More traditional two-parameter sinusoidal models, e.g., PRC(θ) = *A*[sin(*B*) − sin(*B* + θ)], provided less accurate fits (*R*^{2} = 0.8123 ± 0.1009) and were not considered further.

#### Synchrony analyses.

To examine correlation-induced synchrony, MCs were injected with six pairs of noisy inputs, with correlation between inputs of a pair ranging from *C*_{in} = 0 to *C*_{in} = 1 in steps of 0.2. Pairs of inputs were presented to each individually recorded cell once in a pseudorandom order. Resulting binarized spike trains were convolved with a square-wave pulse of 8 ms to account for spike time jitter and differences in spike thresholds (Galán et al. 2008). Convolved spike trains evoked by each input of a correlated pair were then compared within the same cell or between different cells to analyze correlation-induced synchrony within homogeneous or heterogeneous MC populations, respectively. Synchrony was measured using magnitude squared coherence, which normalizes the cross-spectral density between two spike trains by the autospectral density of each individual spike train, yielding a measure of spectral power dependent only on synchrony between spike trains. Cross power spectra were additionally calculated to help visualize synchronization of oscillations. Coherence and power spectra were calculated using the Welch method with a total signal length of 50,000 samples per trace (i.e., 5 s at a sampling rate of 10 kHz), a window size of 1,024 samples, and 50% window overlap. These windowing parameters facilitated clear comparison of output synchrony across conditions while providing sufficient spectral resolution of the target envelope (15–75 Hz). Final results were invariant to square-wave width used in spike train convolution and window size used in spectral analyses (see Fig. 7; Galán et al. 2008). Statistical comparisons made are described in results. For direct comparison with theoretical predictions, phase difference densities were derived from in vitro spike trains by linearly interpolating MC phase across each ISI.

### Simulations

Oscillating MCs recorded in vitro were modeled as simple noisy phase θ oscillators by dθ_{j}/d*t* = ω_{j} + Δ_{j}(θ)[*I*(*t*) + η], where ω_{j} = 2π/*T*_{j} is the natural frequency of neuron *j* calculated using the experimentally derived mean ISI (*T*_{j}), and Δ_{j} is the experimentally derived PRC of neuron *j* (Galán et al. 2008). θ_{j} ranges from 0 to 2π for each model neuron *j*; we take θ_{j} = 0 to be the time of a spike. Model neurons were driven by shared synaptic input, *I*(*t*), identical to in vitro neurons (see above) except that the amplitude was uniformly reduced and balanced with independent white noise η to evoke equivalent “noisy periodic” firing rates and ISI distributions as observed in vitro. As for in vitro neurons, *C*_{in} characterizes the degree of correlation between pairs of frozen colored noise inputs, *I*(*t*). Coherence of output spike trains of model neurons was analyzed as for in vitro neurons (see above). Independent PRC and firing rate manipulations were performed as described in results.

### Mathematical Theory

In this section, we extend previous theory that we have developed to study correlation-induced synchronization of uncoupled neural oscillators. We start with a set of phase models, similar to those above, which arise when the inputs to the oscillators are small. The general model takes the form
*j* = 1, 2, representing two neural oscillators, and ξ_{j}(*t*) is the broadband input of magnitude σ given to the oscillator. θ_{j}, ω_{j}, and Δ_{j} are as defined above. In our prior work (Marella and Ermentrout 2008), we derived equations for the phase difference, ϕ = θ_{2} − θ_{1}, between two identical oscillators (ω_{1} = ω_{2} and Δ_{1} = Δ_{2}) that received a mixture of correlated and uncorrelated white noise *c* is the correlation, η_{1,2} are independent noise sources, and ξ_{C}(*t*) is the common noise]. The result of this calculation is an explicit formula for the density of the phase difference, *P*_{s}(ϕ). This density function can be directly related to the more familiar quantity of spike time cross-correlation, *C*_{ij}(τ), by
*P*_{s}(ϕ) = 1/(2π) is uniform and the cross-correlation is zero. We now report that this formula for the phase difference density can be generalized to take into account other kinds of heterogeneity besides differences in sources of noise, such as heterogeneity in PRC dynamics. Given weak noise (small σ) and small heterogeneities, we obtain the following expression for *P*_{s}(ϕ):
_{j} is the average of the squared amplitude of the PRC, and the factor *ch*(ϕ) provides a sense of how much correlated noise synchronizes the two oscillators (for details, see appendix). Heterogeneity in the natural frequencies has a large effect on the ability of the noisy oscillators to synchronize. If σ^{2} = β (i.e., there is a modest degree of heterogeneity in the frequencies of the 2 oscillators), then *P*_{s}(ϕ) ≈ 1/(2π) − (σ^{2}/2)*ch*′(ϕ)/(2πβ), and *P*_{s}(ϕ) is nearly flat. If β = 0 (i.e., no heterogeneity in the frequencies of the 2 oscillators), then *P*_{s}(ϕ) = *N*/[α_{1} + α_{2} − 2*ch*(ϕ)], where *N* is a normalization constant so that the integral of *P*_{s}(ϕ) is 1. A closed-form solution for *Eq. 5* is possible but not terribly useful, involving a number of integrals that must be evaluated numerically. We thus instead solve *Eq. 5* numerically, using the periodic boundary condition, *P*_{s}(ϕ) = *P*_{s}(ϕ + 2π), and the normalization requirement, ∫_{0}^{2π}*P*_{s}(ϕ)dϕ = 1.

As an illustration of this method applied to experimental data (where oscillators are indeed noisy and heterogeneous), we modeled two oscillators (*Eq. 3*) with ω_{1,2} = 1, using the function form of the PRC (*Eq. 2*) with values derived from two randomly chosen in vitro MCs oscillating at roughly the same periodic rate:

We then applied identical noise to both oscillators (*c* = 1) with σ = 0.25. We integrated the phase models using the Euler method with a step size of 0.05, throwing out the first 10,000 time units and computing 390,000 subsequent time units. We constructed a histogram of ϕ = θ_{2} − θ_{1} between (−π, π) with 100 bins. We also created a “spike train”, *s*_{j}(*t*) = exp(−400{1 − cos[θ_{j}(*t*)]})/0.125371, scaled to be compatible with the binning of ϕ into 100 bins, and computed the cross-correlation of *s*_{1}, *s*_{2}. Finally, we solved *Eq. 5* with β = 0 and *c* = 1 to get *P*_{s}(ϕ). For these noisy oscillators with distinct PRCs (Fig. 1*A*), we note that the numerical phase difference density and spike time correlogram are nearly indistinguishable and that the theoretical density predicted from *Eq. 5* also provides a close match (Fig. 1*B*). This result is consistent across all other oscillator pairs examined. We cannot expect a perfect fit by the theory as the equations are strictly valid only as σ → 0. Nevertheless, the theoretical phase difference density still provides an excellent prediction, thereby relating the phase difference density from dynamical theory to spike time cross-correlation. Moreover, this example demonstrates that we have successfully generalized our previous work (Marella and Ermentrout 2008) to the completely heterogeneous case. We additionally note that, for equal firing rates, the noisy oscillators achieve a large degree of synchrony (Fig. 1*B*) despite distinct PRC dynamics (Fig. 1*A*), although with a noticeable phase shift from ϕ = 0.

Next, we fixed σ = 0.25 while varying β, the difference in natural frequencies, and examined how the shape of the invariant density changed. To characterize the degree of resulting synchrony, we introduce an order parameter: OP := *S* = ∫_{0}^{2π}sin(*x*)*P*(*x*)d*x* and *C* = ∫_{0}^{2π}cos(*x*)*P*(*x*)d*x*. When OP = 0, the density is flat, and when OP = 1, the density has a sharp peak. Thus OP characterizes the degree of synchrony. As stated above, for perfect input correlation (*c* = 1), intrinsic heterogeneity in the PRC dynamics of the two sample oscillators results in a phase shift in the peak of the density *P*_{s}(ϕ). In this example, small increases in β marginally flattened *P*_{s}(ϕ) but surprisingly compensated for the phase shift from ϕ = 0 (Fig. 1*C*). It is thus possible for combined firing rate and PRC heterogeneities to yield a greater proportion of synchronized spikes and thereby encode stimulus-specific information (Friedrich et al. 2004). Morever, this example demonstrates that firing rate and PRC heterogeneities do not exclusively yield reduced synchronization and further suggests that stochastic synchronization of neural oscillators may tolerate a degree of heterogeneity in oscillatory dynamics. In general, however, increasing values of β flatten *P*_{s}(ϕ) and thus reduce OP across all levels of input correlation (Fig. 1*D*).

## RESULTS

### Theoretical Analysis of PRC Heterogeneity

On the basis of the novel theory described above, we predicted the influence of PRC heterogeneity on stochastic synchronization of neural oscillators by using PRCs measured from two in vitro MCs (Fig. 1; see materials and methods). The theory predicted with good accuracy the reduction in synchrony (as measured by the density of phase differences) observed for a given difference in PRCs and intrinsic firing rates. Specifically, we found that as firing rates of two heterogeneous neural oscillators become increasingly different, the phase difference density flattens (Fig. 1, *C* and *D*), corresponding to a decorrelation of the neurons' firing phases. Moreover, as intrinsic firing rates diverge, the peak of the phase difference density shifts (Fig. 1*C*), indicating that pairs of heterogeneous neurons with different intrinsic firing rates will fire with a reliable temporal offset or lag. For any given difference in firing rates, reducing the level of input correlation leads to an approximately linear reduction of output correlation (Fig. 1*D*).

### In Vitro MC PRC Heterogeneity

To investigate the impact intrinsic biophysical diversity has on synchronization of neural oscillations in vitro, we recorded membrane voltages from individual MCs in coronal slices of mouse main OB by using whole cell current-clamp electrophysiology. Each MC was driven by a step current (150–700 pA) added to low-amplitude (σ = 15–20 pA) colored noise to evoke noisy periodic firing (CV_{ISI} < 0.4; Fig. 2*A*) across a range of physiologically relevant rates (15–75 Hz). Critically, this low-CV_{ISI} regime provides a good approximation of periods during which odors evoke periodic firing in MCs and the homologous projection neurons (PNs) of the invertebrate antennal lobe (AL) in vivo (for example, see Cang and Isaacson 2003; Friedrich et al. 2004; Ito et al. 2009; Kashiwadani et al. 1999; Laurent et al. 1996; Margrie and Schaefer 2003; although see discussion). For each oscillating MC, we calculated a STA over the natural period, defined as the cell's mean ISI (Fig. 2*B*). Oscillating MCs exhibited a broad range of STAs (Fig. 2*C*), similar to the diversity measured previously in less periodically firing MCs (Padmanabhan and Urban 2010).

The STA of a neuron firing in a roughly periodic manner (CV_{ISI} < 0.4) can be used to estimate the PRC of the neuron (Ermentrout et al. 2007; Torben-Nielsen et al. 2010). Thus, to translate the biophysical diversity captured by the STA heterogeneity into terms consistent with analyses of oscillatory synchrony, we estimated PRCs for all recorded neurons from their STAs (Fig. 2*B*). Oscillating MCs displayed type II PRC dynamics (Fig. 2*D*) consistent with our previous findings (Ermentrout et al. 2007; Galán et al. 2005). The predominance of type II PRC dynamics in MCs, along with our previous demonstration that stochastic synchronization is enhanced in neurons with type II PRCs (Abouzeid and Ermentrout 2009; Galán et al. 2007; Marella and Ermentrout 2008) agrees well with the observation that odorant stimulation can evoke synchronized oscillatory activity in the rodent OB (Kay et al. 2009). Within the class of type II PRC dynamics, however, MCs exhibited considerable diversity in PRC shape (Fig. 2*D*). To classify in vitro MC PRC heterogeneity, we fit a generalized three-parameter model (*Eq. 2*) to each of the estimated PRCs (Fig. 2*B*; see materials and methods). The proposed PRC model provided an excellent fit to PRCs recorded in vitro (*R*^{2} > 0.99 for all cells). Projection of parameterized PRCs as single points into three-dimensional parameter space (Fig. 2*E*) illustrates the intrinsic diversity and density of MC PRCs recorded in vitro.

Examination of 43 in vitro MC PRCs demonstrated considerable diversity in all model parameters (Fig. 3, *A–C*). Moreover, the parameter values obtained from the fits of these PRCs exhibited low interdependence (Fig. 3, *D–F*). Specifically, the distribution of exponential component parameter values displayed a significant but weak relationship with the distribution of amplitude (Fig. 3*E*; *P* = 5.5 × 10^{−5}, linear regression; *R*^{2} = 0.33) and sinusoidal (Fig. 3*F*; *P* = 0.0051, linear regression; *R*^{2} = 0.18) component parameter values, whereas no significant relationship was observed between amplitude and sinusoidal components (Fig. 3*D*; *P* = 0.62, linear regression). This low parameter interdependence suggests that the proposed three-parameter PRC model is not unnecessarily complex and further supports the conclusion that PRC heterogeneity is an intrinsic property of in vitro MCs that cannot be attributed to variation in any single parameter.

### Firing Rate-Dependent Modulation of PRC Dynamics

Changes in periodic firing rates can modulate PRC type and PRC shape within a single dynamical type (Fink et al. 2011; Gutkin et al. 2005; Schultheiss et al. 2010; Stiefel et al. 2008, 2009). Thus we next asked, to what extent is the observed MC PRC diversity a product of firing rate-dependent modulation? To address this question, we first characterized the range of noisy periodic firing rates evoked in the population of recorded MCs by step current injections overlaid with low-amplitude colored noise. For any given magnitude of step current, the population of MCs exhibited a broad range of resulting firing rates (Fig. 4*A*). Indeed, no significant relationship was found between step current magnitude and evoked noisy periodic firing rate across the MC population (*P* = 0.18, linear regression). This heterogeneity in evoked firing rates was consistent with the substantial diversity observed within MC FI relationships (Fig. 4*B*; Padmanabhan and Urban 2010), and could not be directly attributed to differences in input resistance (*P* = 0.53, linear regression). In vitro MC heterogeneity is thus evident not only in the cells' PRCs but also in differences in preferred firing rate ranges.

To assess how the heterogeneity in evoked noisy periodic firing rates contributes to the heterogeneity in calculated PRCs, we next quantified the difference in parameterized PRC components between MCs firing at different periodic rates. Specifically, we calculated the mean absolute difference in PRC model amplitude, sinusoidal node, and exponential component parameters between MCs oscillating with firing rate differences of 0–5, 5–10, 10–15, and 15–20 Hz (Fig. 4, *C–E*). Across the MC population, only the exponential component of parameterized PRCs exhibited any significant firing rate dependence (Fig. 4*E*; *P* = 0.0021, 1-way ANOVA), with MCs oscillating at highly different rates (15–20 Hz) exhibiting larger differences in exponential components than MCs oscillating at highly similar rates (0–5 Hz). Differences in individual parameterized PRC components are thus largely independent of firing rate differences between MCs.

Next, we visualized full PRC models and clustering of PRCs in three-dimensional parameter space on the basis of firing rate (Fig. 4, *F1* and *F2*). Consistent with the above analysis, the set of parameterized PRCs exhibited a general trend toward greater phase delays (Fig. 4*F1*) and lower exponential component parameters (Fig. 4*F2*) with higher firing rates. Note, however, that there still exists considerable heterogeneity in PRCs recorded from MCs firing at identical rates. To provide an additional objective measure for how PRC and firing rate heterogeneities relate, we performed multiple linear regression between parameterized PRCs and recorded periodic firing rates (Fig. 4*F3*). A significant (*P* = 0.0038, 1-way ANOVA) but weak (*R*^{2} = 0.29) relationship was found between the recorded firing rate and the rate predicted by a linear combination of PRC parameters, with only the exponential component providing significant predictive power (*P* = 0.0014, post hoc Tukey's test). In other words, 29% of the heterogeneity in recorded PRCs could be attributed to the heterogeneity in evoked periodic firing rates. Thus PRC and firing rate heterogeneity are largely independent across the population of recorded MCs.

We additionally considered to what extent variance in spiking periodicity and passive membrane properties related to the recorded PRC heterogeneity. As expected, greater periodicity (i.e., lower CV_{ISI}) was associated with lower amplitude PRCs (Fig. 4, *G1* and *G2*). In other words, MCs exhibiting reduced responsiveness to the aperiodic noisy current injection demonstrated behavior closer to that of a pure oscillator (CV_{ISI} = 0). A general trend toward lower exponential component parameters with greater periodicity was also observed (Fig. 4*G2*), possibly due to the weak trend toward lower CV_{ISI} with higher firing rates (data not shown; Padmanabhan and Urban 2010). In total, heterogeneity across the set of parameterized PRCs was strongly (*R*^{2} = 0.52) and significantly (*P* = 2.2 × 10^{−6}, 1-way ANOVA) related to the observed variance in recorded CV_{ISI} (Fig. 4*G3*), with both amplitude (*P* = 0.020, post hoc Tukey's test) and exponential (*P* = 0.0074, post hoc Tukey's test) component parameters providing significant predictive power. In contrast to heterogeneity in periodic firing rates, PRC heterogeneity exhibited a strong (*R*^{2} = 0.56) and significant (*P* = 3.6 × 10^{−7}, 1-way ANOVA) relationship to variance in recorded input resistance across the MC population (data not shown), as expected from their explicit relationship (*Eq. 1*). Variance in recorded membrane time constant was only weakly related to the observed heterogeneity in PRCs (*P* = 0.036, 1-way ANOVA; *R*^{2} = 0.19), whereas no relationship existed between PRC heterogeneity and the magnitude of step current injection used to evoke periodic firing (Fig. 4, *H1–H3*; *P* = 0.12, 1-way ANOVA), consistent with intrinsic biophysical properties underlying PRC (and firing rate) heterogeneity.

To explore PRC modulation by firing rate changes further, we drove a subset of MCs at multiple firing rates and quantified the sensitivity of parameterized PRC within each MC across different firing rates. Using this within-cell approach, we observed a shift toward lower exponential component parameters with increased firing rates (Fig. 5*A*), congruent with our analysis across the MC population (Fig. 4*F2*). In addition, we observed lower PRC amplitudes with increased firing rates (Fig. 5*A*). Quantification confirmed that both PRC model amplitude and exponential parameters significantly decreased with increased firing rate (Fig. 5*B*). Collectively, these firing rate dependencies increased type II PRC dynamics by reducing positive (phase advance) PRC values and further decreasing negative (phase delay) PRC values (Fig. 5*C*). The net effect of firing rate-dependent modulation of MC PRCs proved relatively small over a large range of firing rates (15–75 Hz), however (Fig. 5*C*), and thus cannot account for the much larger degree of total MC PRC heterogeneity. Intrinsic biophysical diversity within in vitro MCs thus gives rise to substantial and largely independent PRC and firing rate heterogeneity, which we collectively refer to as heterogeneity in “oscillatory dynamics.”

### Effect of Heterogeneity in Oscillatory Dynamics on Correlation-Induced Synchronization In Vitro

Our data thus far establish a considerable degree of intrinsic biophysical diversity within oscillating MCs, as measured by STA and PRC heterogeneity and reflected in heterogeneous FI relationships. We next sought to investigate the functional consequences of this cellular diversity. In particular, how does intrinsic biophysical diversity affect correlation-induced synchronization? To address this question, we drove a subset of 27 MCs with pairs of low-amplitude noisy correlated input (Fig. 6), mimicking varying degrees of shared glomerular and granule cell synaptic input, and calculated the resulting correlation-induced coherence between evoked spike trains as a measure of synchronization. To eliminate possible contributions of coupling-induced synchronization via lateral synaptic inhibition to our measures of synchrony, we recorded from individual MCs in different slices. This in vitro approach, similar to that used by Galán et al. (2006) and Markowitz et al. (2008), allowed us to then simulate uncoupled pairs of real MCs receiving correlated input and to remove the potential influence of local circuitry on correlation. To specifically isolate the effect of intrinsic MC diversity on the resulting correlation-induced synchrony, we first measured the coherence between spike trains evoked within the same MC by correlated input (thus simulating a “homogeneous” pair of MCs) and then compared this with the coherence measured between spike trains evoked in two different MCs by correlated input (thus simulating a “heterogeneous” pair of MCs) (Fig. 6). Importantly, the weak fluctuations used (σ = 15 pA) do not reliably drive spike timing (Galán et al. 2008), allowing us to examine spike correlation in the oscillatory regime and differentiating our experiments from those focused specifically on noise-induced spike precision and reliability (Mainen and Sejnowski 1995). Critically, differences in periodic firing rates can significantly affect stochastic synchrony (Markowitz et al. 2008; and see below). Thus, to distinguish effects of biophysical diversity (as measured by PRC heterogeneity) from effects of firing rate differences, we first focused on correlation-induced spike time coherence between MCs firing at highly similar (|ΔFR_{ij}| ≤ 5 Hz) periodic rates.

Foremost, we observed a clear effect of input correlation on spike time cross power spectra and coherence in both homogeneous (Fig. 7, *A* and *B*) and heterogeneous (Fig. 7, *C* and *D*) MC populations, recapitulating the phenomenon of correlation-induced synchronization (Galán et al. 2006). This was further evident in mean spike time correlogram peaks, once time was rescaled to the natural period of oscillations (Fig. 7, *A* and *C*, *insets*). To quantitatively compare the level of correlation-induced synchronization in homogenous vs. heterogeneous (|ΔFR_{ij}| ≤ 5 Hz) MC populations, we integrated the spike time coherence from 15–75 Hz for each level of input correlation considered. With the use of this measure, the extent of synchronization was significantly greater in homogeneous than in heterogeneous MC populations (Fig. 7*E1*), with cell heterogeneity imposing up to a 30% reduction in maximum output coherence (Fig. 7*E2*). Intrinsic biophysical diversity within MCs oscillating at similar rates thus significantly limits (but does not abolish) stochastic synchronization.

Markowitz et al. (2008) previously demonstrated that small firing rate differences can have a strong effect on output synchrony of pyramidal neuron pairs driven by perfectly correlated input. To what extent this firing rate effect extends across a population of neurons, depends on the level of input correlation, or relates to physiological neuronal heterogeneity remains unknown, however. Thus we next quantified the effect of firing rate divergence on stochastic synchrony. To do so, we compared the spike time coherence between oscillating MCs with firing rate differences of 0–5, 5–10, 10–15, or 15–20 Hz throughout the range of physiological beta/gamma frequencies. As expected, increasing firing rate divergence significantly reduced the extent of stochastic synchronization observed (Fig. 7*F1*), with each 10-Hz increase in firing rate difference imposing up to a 25% reduction in maximum output coherence (Fig. 7*F2*). This effect could not be attributed to any systematic shift in biophysical diversity, since PRC differences were consistent across all MC comparisons (Fig. 7*G*). Moreover, firing rate differences within homogeneous MC populations also significantly reduced the extent of correlation-induced synchronization (Fig. 7*H*) despite low firing rate-dependent PRC modulation (Fig. 5, *B* and *C*). Thus heterogeneity in oscillatory dynamics, comprising both PRC and firing rate heterogeneity, plays a key role in governing stochastic synchronization. Specifically, physiological PRC heterogeneity and low to moderate firing rate differences each separately reduce the maximum synchrony possible by up to 25–30%. Combined, these two sources of heterogeneity impose up to a 40% reduction in maximum output synchrony, consistent with a small but significant fraction of PRC heterogeneity (and its impact on output synchrony) arising from firing rate differences (see above). We additionally note that, surprisingly, biophysically diverse MCs oscillating at rates as different as 20 Hz still supported a measure of correlation-induced synchronization (Fig. 7*F1*). Our results thus establish stochastic synchronization as a robust phenomenon of real neurons in vitro.

### Isolating the Effect of PRC Heterogeneity and Firing Rate Differences on Correlation-Induced Synchronization

We next sought to determine how completely PRC heterogeneity accounted for the reduction in stochastic synchrony observed between neurons oscillating at highly similar rates. That is, we sought to test whether aspects of neuronal dynamics not captured by the PRC contributed significantly to the reduced synchrony observed between vs. within MCs. We thus simulated MC responses using simple models in which all of the neuronal dynamics are derived from the PRC. Specifically, for each in vitro MC, we constructed a noisy phase oscillator using the MC's recorded oscillatory frequency and parameterized PRC (see materials and methods). We then performed analyses of synchrony equivalent to those performed on the spike trains from the real neurons above. Within this highly reduced framework, we observed levels of spike time cross power spectra and correlation-induced spike time coherence in homogeneous (Fig. 8, *A* and *B*) and heterogeneous (Fig. 8, *C* and *D*) oscillator populations comparable to those in the corresponding in vitro populations (Fig. 7, *A–D*). Critically, the PRC heterogeneity measured in vitro proved sufficient in these simulations to account for at least the reduction in synchrony observed among heterogeneous MCs in vitro (compare Fig. 8*E* with Fig. 7*E1*). In other words, differences in firing rates were not required to explain the experimentally observed reduction in synchrony. Likewise, small to moderate differences in oscillatory frequencies yielded at least the same reduction in maximum output coherence among noisy phase oscillators (Fig. 8*F1*) as among in vitro MCs (Fig. 7*F*). Our simulations thus confirmed that the range of PRC heterogeneity and firing rate differences observed in vitro alone are sufficient to account for the observed changes in stochastic synchronization. Reducing PRC differences to zero significantly enhanced synchronization within heterogeneous populations of simulated oscillators with low firing rate differences (0–5 Hz) to levels near those observed within homogeneous populations (Fig. 8*E*). In heterogeneous populations with high firing rate differences (15–20 Hz), however, reducing PRC differences to zero had no significant effect on synchronization of the simulated oscillators (Fig. 8*G*). Our simulations thus further revealed that *1*) both PRC heterogeneity and firing rate differences limit stochastic synchronization given low to moderate firing rate differences, and *2*) high firing rate differences have a dominant effect on stochastic synchronization.

The effect of biophysical diversity on stochastic synchrony is illustrated in an example in which we examined the response of a random subset of in vitro MCs with heterogeneous PRCs (Fig. 9*A*) firing at highly similar rates to a series of highly correlated inputs (*C*_{in} = 0.8). Spike time raster plots (Fig. 9*B*) show that each MC exhibited clear noisy periodic activity on each sweep as well as a considerable degree of sweep-to-sweep synchrony induced by the small-amplitude, highly correlated input (Fig. 9*C*). Moreover, the subset of MCs also exhibited clear epochs of between-cell synchrony, often with different MC pairs aligning at different times, producing a pattern of phasic synchronization and yielding greater within-cell than between-cell oscillatory synchrony across the entire MC subset at any given instance. This effect was particularly evident in peristimulus time histograms (Fig. 9*D*, *top*), where firing probability within a single cell exhibited sharper peaks and deeper interspike valleys than the firing probability across the entire subset (reflected in the absolute probability difference calculated in Fig. 9*D*, *bottom*). Similar investigation in our modeling framework with the use of a large population of noisy phase oscillators yielded equivalent results (Fig. 9, *E–H*).

### Impact of Single PRC Component Variation on Correlation-Induced Synchronization

Intrinsic biophysical diversity within MCs in vitro, as measured by PRC heterogeneity, thus significantly affects correlation-induced synchronization of neural oscillations (Figs. 7⇑–9). This heterogeneity within MC PRCs is captured in our parameterized PRC model and arises from variance in three components: overall amplitude, sinusoidal node, and the balance between phase advance and delay regimes (Fig. 3). In our next experiment, we thus examined how diversity within each of these PRC components independently impacts correlation-induced synchronization. Understanding which aspects of PRC shape have the greatest effect on synchronization will ultimately help identify the biophysical properties underlying the overall effect of cellular diversity on synchronization.

We simulated a homogeneous population of noisy phase oscillators with firing rates and model PRC parameters matching the mean values observed in vitro. We then varied a single PRC parameter throughout its observed in vitro range (Fig. 10, *A*, *D*, and *G*), calculated the net difference in resulting PRCs between oscillators (Fig. 10, *B*, *E*, and *H*), and integrated the resulting spike time coherence from 15–75 Hz between each pair of simulated oscillators for perfectly correlated input (Fig. 10, *C*, *F*, and *I*). It is important to note that this simple approach neglects the significant interdependence observed among the distribution of PRC component parameter values for the population of recorded MCs (Fig. 3). Given the relatively weak (*R*^{2} ≤ 0.33) relationships observed among these parameters, however, results of this analysis should still largely hold for MCs and generalize to any population of neurons exhibiting noisy periodic activity (CV_{ISI} < 0.4) and largely or wholly independent variance among the three PRC components.

In keeping with our in vitro results, nearly all resulting PRCs were clearly type II, with the exception of the highest values of either the sinusoidal (Fig. 10*D*) or exponential (Fig. 10*G*) component parameters. In general, the greatest degree of synchrony was observed along the diagonal of each coherence plot (Fig. 10, *C*, *F*, and *I*), which corresponds to the dark blue diagonal demarcating zero PRC difference in Fig. 10, *B*, *E*, and *H*. This reaffirms the impact of PRC heterogeneity on correlation-induced synchronization. However, neither the degree of synchrony observed between identical oscillators nor the effect of PRC differences on synchrony was uniform across the range of parameter values considered. For variation in the amplitude component parameter, small levels of PRC heterogeneity strongly desynchronized oscillators with low-amplitude PRCs (Fig. 10*C*, *bottom left* corner), whereas synchronization of oscillators with high-amplitude PRCs tolerated substantially greater PRC heterogeneity (Fig. 10*C*, *top right* quadrant). For variation in the sinusoidal component parameter, high levels of synchrony were only observed between oscillators with similarly low sinusoidal component parameters (Fig. 10*F*, *bottom left* quadrant), even among oscillators with highly similar PRCs (Fig. 10*F*, *top right* corner) and despite a reduced effect of parameter variation on total PRC differences at high sinusoidal component parameters (Fig. 10*E*). Because high values of the sinusoidal component parameter shift the balance between phase advance and delay (Fig. 10*D*), this result supports the facilitation of synchronization by type II vs. type I PRC dynamics (Abouzeid and Ermentrout 2009; Galán et al. 2007; Marella and Ermentrout 2008). Variation in the exponential component parameter, which also determines the balance between phase advance and delay (Fig. 10*G*), had a surprisingly smaller effect on stochastic synchronization (Fig. 10*I*). Tolerance of synchrony to variation in the exponential component parameter likely arose due to *1*) the minimal impact parameter variation had on total PRC differences (Fig. 10*H*) and *2*) the reduction in overall PRC amplitude and window of responsive phases (Fig. 10*G*). Altogether, the impact of PRC heterogeneity on stochastic synchronization thus most strongly depends on PRC type and amplitude.

### Comparison of In Vitro and Simulation Results With Theoretical Predictions

In keeping with our theoretical predictions (see materials and methods), our in vitro and simulation results thus far demonstrate that both PRC heterogeneity and firing rate divergence significantly impact correlation-induced synchronization of neural oscillators. To provide more direct comparison between our experimental results and theoretical predictions, we also examined phase difference densities and spike time correlograms from our in vitro and simulation data (Fig. 11). Instantaneous MC phase was estimated from in vitro voltage traces by linearly interpolating phase across each ISI (Fig. 11*A*). Shown clearly for two random in vitro MCs with heterogeneous PRCs (Fig. 11*B*) and similar firing rates given perfectly correlated input, the in vitro phase difference density and spike time correlogram closely agreed (Fig. 11*C*), revealing a peak in the distribution with a slight phase shift from ϕ = 0 indicative of considerable oscillatory synchrony similar to our theoretical predictions for two other random sample MCs (Fig. 1*B*). Furthermore, the same result was derived by solving our novel analytic theory (*Eq. 5*) for the phase difference density (Fig. 11*C*) using model PRCs parameterized from the sample MCs (Fig. 11*B*). This example demonstrates that our theoretical framework can accurately capture the impact biophysical diversity among real neurons has on synchronization of periodic firing.

Given this agreement between our mathematical theory and experimental observations, we would expect that the reduction in synchrony imposed by cellular heterogeneity and firing rate differences among in vitro MCs (Fig. 7) should also manifest in a flattening of the in vitro phase difference density, as predicted by our theory (Fig. 1). Across the population of 27 MCs for which analyses of oscillatory synchrony were performed, we first observed that the reduction in maximum output coherence imposed by intrinsic biophysical diversity (Fig. 7*E2*) indeed translated into a measurable flattening of the average phase difference density and spike time correlograms (Fig. 11*D*). Moreover, we again note the close qualitative relationship observed between the phase difference density and spike time correlogram (Fig. 11*D*), as established by our theory (*Eq. 4*). Our estimate of in vitro MC phase also recapitulated the phenomenon of correlation-induced synchrony in terms of phase difference densities, shown in Fig. 11*E* for the homogeneous MC population, as well as the pronounced (although not complete) flattening of the phase difference density with increasing firing rate divergence (Fig. 11*F*), as predicted by our theory (Fig. 1, *C* and *D*). Equivalent results were also found using our reduced neural simulations (Fig. 11, *G–I*). In total, we thus observed close agreement among our theoretical predictions, which are in terms of phase difference densities, and our electrophysiological and simulation data, lending considerable credance to our core finding that both PRC and firing rate heterogeneity place important constraints on stochastic synchronization.

## DISCUSSION

### Summary

Analysis of the mechanisms of oscillatory synchrony has been an important contribution of mathematical theory and modeling to neuroscience. Although the impact of heterogeneity on coupling-induced oscillatory synchrony has been extensively studied, models of correlation-induced synchrony (including our own previous work) have typically not taken into account the fact that neurons, even those of a single molecular type, are not identical to each other but rather differ in their intrinsic properties (Marder and Goaillard 2006). Here, for the first time, we examine the degree to which biophysical diversity and firing rate differences in a population of real neurons influence correlation-induced oscillatory synchrony. The neuronal heterogeneity examined, which may be due to differences in ion channel expression, membrane properties, and morphology, will also impact the neuronal response to stimuli and the coding properties of populations (Padmanabhan and Urban 2010). Here, we first measured the degree of diversity in a specific population of neurons (main OB MCs) recorded in vitro via whole cell current-clamp electrophysiology. We then analyzed the output of these neurons to determine how the measured biophysical diversity influences oscillatory synchronization of MCs in response to correlated fluctuating inputs, such as might be received from common synaptic partners. Specifically, we used MC recordings to estimate PRCs and thus quantified the response of periodically firing neurons to transient inputs. We found that MCs exhibit a broad range of type II PRCs that could be well characterized by a relatively simple phenomenological equation with three parameters. From our in vitro data, we estimated the density of MCs in the space defined by these PRC parameters and used this estimate to predict how MC diversity impacts correlation-induced oscillatory synchrony. Collectively, our experimental findings, together with the results of simulations based closely on the recordings of in vitro MCs, provide good agreement with predictions of novel analytic theory also described here.

### Heterogeneity in Mechanisms of Oscillatory Synchrony

Previous work has demonstrated that heterogeneity in networks of oscillating neurons can reduce synchrony induced by synaptic or gap junction coupling (Golomb and Rinzel 1993; Wang and Buzsáki 1996; for review, see Kopell and Ermentrout 2002). Coupling in such networks has generally been modeled as “all-to-all” such that each cell is coupled to every other cell. In a majority of cases, heterogeneity reduced the ability of the neurons to synchronize, with large networks often undergoing a phase transition to a synchronized state as heterogeneity is reduced (Kuramoto 2003). Like our present results, the way in which neurons synchronize (such as via synaptic inhibition) can make the neuronal synchrony more or less susceptible to heterogeneity (Chow 1998; White et al. 1998). Gap junctions, through their ability to average the potentials of connected neurons, can reduce the effects of heterogeneity (Kopell and Ermentrout 2004). Even when neurons have identical biophysical properties, heterogeneity in the number of inputs can impair oscillatory synchrony (Skinner et al. 2005). Analysis of the effects of heterogeneity on correlation-induced synchronization of oscillating neurons has previously been limited to study of the effects of firing rate and input heterogeneity on transient synchronization of neurons (Brody and Hopfield 2003; Markowitz et al. 2008).

### Heterogeneity of Mitral Cells as Oscillators

Our previous work has shown that MCs are heterogeneous in their intrinsic properties and that this source of cellular diversity significantly impacts encoding of stimulus features (Padmanabhan and Urban 2010). Physiological heterogeneity will also impact the ability of real neurons to synchronize. To analyze these effects, we have treated neurons (in this case, MCs) as noisy oscillators and considered how two different physiological components of heterogeneity impact correlation-induced neural synchrony. We examined heterogeneity in intrinsic firing rate or natural frequency and heterogeneity in the neuronal response to input, as captured by the STA and PRC. Because we specifically examined oscillatory synchrony among periodically firing neurons, the diversity in neuronal response to input was best captured by differences in PRCs. PRC shape and magnitude in homogeneous populations significantly influence correlation-induced synchronization (Abouzeid and Ermentrout 2009; Galán et al. 2007; Marella and Ermentrout 2008). We now demonstrate that both firing rate and PRC heterogeneity place important constraints on stochastic synchronization. Specifically, physiological diversity encompassed within PRC heterogeneity resulted in up to a 30% reduction in the maximum output synchrony in vitro, whereas moderate firing rate differences (∼10 Hz) similarly imposed up to a 25% reduction. Combined, these two sources of heterogeneity yielded up to a 40% reduction in the maximum output synchrony in vitro. There are a number of ways to understand the impact of heterogeneity in this context. From the standpoint of generating synchronous spiking in a heterogeneous population, up to a ∼30% larger population would be required to generate the same number of synchronized spikes as a homogeneous population, and an even larger population would be required given moderate firing rate differences among the heterogeneous neurons. Alternatively, the total impact of intrinsic biophysical diversity among this class of principal neurons is approximately equivalent to a decrease in input correlation by as much as 0.2. This is similar to the range of input correlations observed between neurons in a number of brain areas (Cohen and Kohn 2011), indicating a large role for cell-to-cell heterogeneity in regulating oscillatory synchrony.

Given this result, we predict that the degree of PRC heterogeneity among MCs (or among any cell type engaged in periodic firing) will significantly impact the degree of stochastic synchronization possible for a given level of firing rate heterogeneity. Such changes to PRC heterogeneity could arise through activity-dependent changes in channel expression, internalization and recycling, and posttranslational modifications or through changes in channel (in)activation kinetics via neuromodulation and could occur largely independently of (or on a slower timescale than) moment-to-moment changes in driving synaptic input and consequent periodic firing rate. Such a mechanism could thereby enhance synchrony among neurons subject to similar activity levels or neuromodulatory input and decrease synchrony among differentially activated or modulated neurons. Modifying neuronal diversity thus provides a novel mechanism of modulating output synchrony without necessarily altering the average properties of the neurons involved. For example, cholinergic input, which modulates numerous voltage-dependent channels in the hippocampus and cortex and strongly facilitates synchronous network oscillations (Cobb and Davies 2005), can also homogenize PRCs into a single dynamical type (Stiefel et al. 2008). Conversely, prolonged febrile seizures in rodents precipitate a long-term increase in the variance (but not mean) of resting membrane potential of CA1 stratum oriens interneurons (Aradi and Soltesz 2002), likely by altering overall input and leak conductances. Such diversification of interneuronal responsiveness should in turn limit the degree of oscillatory synchrony possible and thus may serve as a homeostatic response to epileptic network rhythms. Modulation of oscillatory synchrony due to altered average properties across a population may thus be enhanced or countered by modifications in the variance of those properties. Such changes in neuronal heterogeneity may even contribute to the pathological alterations in oscillatory synchrony observed in a variety of neurological disorders, including schizophrenia, epilepsy, and autism (Uhlhaas and Singer 2006).

### Biophysical Components of PRC Heterogeneity

In this study, we employed a three-parameter phenomenological model to classify in vitro MC PRC heterogeneity. This simple model affords multiple advantages in the study of how cellular diversity impacts correlation-induced synchrony. First, the proposed model makes no assumptions about underlying ion channel distributions and thus can be generalized to any type of neuron firing in a roughly periodic manner. Second, the proposed model can reproduce the periodicity displayed by many neuronal PRCs while minimizing the potential of overfitting phasic perturbation data with an expanded Fourier series (Galán et al. 2005; Netoff et al. 2012; Torben-Nielsen et al. 2010). Third, the proposed model enables PRC heterogeneity to be decomposed into variation in tractable PRC features (i.e., amplitude, sinusoidal node, and phase delay/advance balance, rather than abstract spline or harmonic coefficients) and then related to other cellular properties, such as firing rate. Finally, heterogeneity in each of these tractable model features can be directly related to differences in conductance and/or time constants of ion channels to provide critical biophysical insight into cellular diversity. For example, increasing the activation rate of delayed rectifying potassium channels or increasing the inactivation rate of sodium channels selectively amplifies phase delays in a firing rate-dependent manner (Fink et al. 2011), equivalent to a decrease in the exponential component parameter of the PRC model at low sinusoidal component parameter values. Likewise, increasing the conductance of slow low-threshold potassium currents can shift the node between phase delay and advance regimes to later phases while selectively reducing the amplitude of phase advances (Ermentrout et al. 2001; Gutkin et al. 2005; Pfeuty et al. 2003; Stiefel et al. 2008, 2009), equivalent to a decrease in the sinusoidal component parameter of the PRC model coupled to an increase in the exponential component parameter. In general, enhancing depolarizing currents yields greater type I PRC dynamics (equivalent to increasing sinusoidal and exponential component parameters), whereas enhancing hyperpolarizing currents yields greater type II PRC dynamics (equivalent to decreasing sinusoidal and exponential component parameters) (for review and further examples, see Ermentrout et al. 2012). Modulating a cell's electrotonic compactness, such as by changing input resistance, provides the simplest (although not only) means of purely modulating PRC amplitude.

### Conserved Contributions of Heterogeneity in Oscillatory Dynamics to Olfactory Encoding

Beyond providing a plastic substrate for influencing neuronal synchrony, how might heterogeneity in oscillatory dynamics contribute to olfactory processing in general? During the initial stages of olfaction, odors trigger widespread synchronized oscillatory activity in the OB of mammals (Adrian 1942; Eeckman and Freeman 1990; Kashiwadani et al. 1999) and lower vertebrates (Friedrich and Laurent 2001), and the homologous AL of invertebrates (Ito et al. 2009; Laurent and Davidowitz 1994; Stopfer et al. 1997; Tanaka et al. 2009). Within this network rhythm, odors can drive transient noisy periodic firing in MC and PN (for example, see Cang and Isaacson 2003; Friedrich et al. 2004; Ito et al. 2009; Kashiwadani et al. 1999; Laurent et al. 1996; Margrie and Schaefer 2003), with a fraction of MC and PN spikes phase-locked to the network local field potential (Buonviso et al. 2003; Friedrich et al. 2004; Ito et al. 2009; Kashiwadani et al. 1999; Laurent and Davidowitz 1994; Tanaka et al. 2009). Each odor evokes a stimulus- and concentration-specific pattern of phasic MC/PN activation and oscillatory synchronization that progressively diverges across time from initial glomerular input patterns (Bathellier et al. 2008; Friedrich et al. 2004; Friedrich and Laurent 2001; Niessing and Friedrich 2010; Stopfer et al. 2003; Wehr and Laurent 1996; Yaksi and Friedrich 2006), thereby encoding stimulus-specific information while the total change in firing rate across the MC/PN population is relatively constant and odor-independent (Kay and Stopfer 2006). Moreover, despite the odor and concentration independence of synchronized MC/PN spike phase within the network local field potential (Friedrich et al. 2004; Laurent and Davidowitz 1994; Stopfer et al. 2003; Wehr and Laurent 1996), the proportion of synchronized spikes each MC fires can be odor dependent and uncorrelated with the total odor-evoked MC firing rate (Friedrich et al. 2004), thereby providing a second channel of stimulus-specific information. Heterogeneity in the propensity of oscillating MCs to synchronize (i.e., PRC diversity) thus may directly contribute to at least two facets of olfactory encoding: first, by diversifying the order with which MCs synchronize their activity (and thereby propagate activity to higher brain regions) to construct odor-specific activity patterns; and second, by modulating the fraction of synchronized spikes each MC fires, as the PRC strictly determines the timing (and not number) of spikes fired. Firing rate-dependent modulation of MC PRC dynamics could also provide an additional variable by which stimulus-specific information may be conveyed to downstream targets via synchronized oscillatory activity. Given that MCs are subject to both correlated input from common feedforward and lateral interactions, as well as inhibitory coupling via granule cells, how coupling- and correlation-induced mechanisms of oscillatory synchrony interact in olfactory encoding will ultimately prove a key area of future investigation. As a caveat, it is important to note that odors can also evoke aperiodic (CV_{ISI} > 0.4) firing in MCs (for example, see Friedrich and Laurent 2004; Kay and Laurent 1999; Rinberg et al. 2006) and that the full role of oscillatory synchrony in olfaction remains uncertain.

## GRANTS

This work was supported by fellowships from the Achievement Rewards for College Scientists Foundation and R. K. Mellon Foundation (to S. D. Burton), National Institute on Drug Abuse Predoctoral Training Grant 2T90 DA022762-06 (to S. D. Burton), and National Institute on Deafness and Other Communication Disorders Grant 5R01 DC011184-06 (to G. B. Ermentrout and N. N. Urban).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## AUTHOR CONTRIBUTIONS

S.D.B., G.B.E., and N.N.U. conception and design of research; S.D.B. and G.B.E. performed experiments; S.D.B., G.B.E., and N.N.U. analyzed data; S.D.B., G.B.E., and N.N.U. interpreted results of experiments; S.D.B., G.B.E., and N.N.U. prepared figures; S.D.B., G.B.E., and N.N.U. drafted manuscript; S.D.B., G.B.E., and N.N.U. edited and revised manuscript; S.D.B., G.B.E., and N.N.U. approved final version of manuscript.

## APPENDIX

#### Relating the Phase Difference Density to Spike Time Correlation

Let *Y*(θ) be a pulse-like function that converts the phase of a neural oscillator to a spike [e.g., *Y*(θ) = δ(θ), where δ(θ) is the Dirac delta function]. We implicitly make the function periodic. Thus we create “spike trains,”*y*_{j}(*t*) = *Y*[θ_{j}(*t*)], from the phases (θ_{j}) of the oscillators. The cross-correlation (with mean firing rate subtracted) is then
*T* → ∞ this is an average, and thus

By definition, θ_{j} − θ_{i} = ϕ so that to lowest order, θ_{j}(*t* + τ) ≈ θ_{j}(*t*) + τ ≈ θ_{i}(*t*) + τ + ϕ (where we assume that the frequency of the oscillator is 1). The probability Pr(θ_{i}, θ_{j}) can be expressed in terms of θ_{i} and ϕ = θ_{j} − θ_{i}. For small noise, the density is uniform in θ_{i} so that

For *Y*(θ), a Dirac delta function, we can explicitly obtain

#### Deriving the Phase Difference Density for Two Heterogeneous Uncoupled Oscillators Driven by Correlated Noise

To formally explore the impact of oscillator heterogeneity on correlation-induced synchronization of uncoupled oscillators, we start with two white noise-driven uncoupled oscillators under the assumptions that the noise and heterogeneities are small. After reduction to a pair of phase equations (Teramae et al. 2009), we obtain the Ito differential equations:
*a*_{j}(θ) is heterogeneity in, say, currents, Δ_{j}(θ_{j}) are the PRCs, *W*(*t*) is the common (correlated) part of the white noise stimulus and *W*_{j}(*t*) is the independent part of the noise, and σ is the magnitude of the partially correlated noise. What is the stationary density of the phase or time difference, *P*_{s}(θ_{2} − θ_{1})? To find this, we make a simple change of coordinates. Let η = θ_{1} and let ϕ = θ_{2} − θ_{1}. We can then rewrite the stochastic equations as

Following Gardiner (1985), we can write the associated Fokker-Planck equation for the probability density, *P*(η, ϕ, *t*). Because we are only interested in the steady state (∂*P*/∂*t* = 0), we obtain the following equation:

We want to obtain a simplified expression for the marginal density, *P*_{s}(ϕ) := ∫_{0}^{2π}*P*(η, ϕ)dη. For weak noise (σ small) and small heterogeneities (*a*_{j} small), the density of η is nearly uniform, so we can approximate the full density as *P*(η, ϕ) = *P*_{s}(ϕ)/(2π). With this approximation, we can integrate the Fokker-Planck equation with respect to η and use the fact that the boundary conditions are periodic to obtain a simple differential equation for the marginal density:
^{2}α, where σ is the magnitude of the noise. We note that α_{1} + α_{2} ≥ 2*h*(ϕ) with equality only if Δ_{1}(θ) = Δ_{2}(θ) and ϕ = 0. This is explained by the following: ∫[*f*(*x*) − *g*(*x*)]^{2}d*x* ≥ 0 with equality only if *f*(*x*) = *g*(*x*). Expanding this yields ∫*f*(*x*)^{2}d*x* + ∫*g*(*x*)^{2}d*x* ≥ 2∫*f*(*x*)*g*(*x*)d*x*. Substituting *f* = Δ_{1}(*x*) and *g* = Δ_{2}(*x* + ϕ) gives ∫Δ_{1}(*x*)^{2}d*x* + ∫Δ_{2}(*x* + ϕ)^{2}d*x* ≥ 2∫Δ_{1}(*x*)Δ_{2}(*x* + ϕ)d*x*.

- Copyright © 2012 the American Physiological Society