Patterns of interval correlations in neural oscillators with adaptation

1. Introduction

The nerve cells of the brain are complex physical systems. They generate action potentials (spikes) by a nonlinear, adaptive, and noisy mechanism. In order to understand signal processing in single neurons, it is vital to analyze the sequence of the interspike intervals (ISIs) between adjacent action potentials. There is experimental evidence accumulating that the spiking in many cases is not a renewal process, i.e., a spike train with mutually independent ISIs, but that intervals are typically correlated over a few lags (Lowen and Teich, 1992; Ratnam and Nelson, 2000; Neiman and Russell, 2001; Nawrot et al., 2007; Engel et al., 2008) [further reports are reviewed in (Farkhooi et al., 2009; Avila-Akerberg and Chacron, 2011)]. These correlations are a basic statistics of any spike train with important implications for information transmission and signal detection in neural systems (Ratnam and Nelson, 2000; Chacron et al., 2001, 2004; Avila-Akerberg and Chacron, 2011) and man-made signal detectors (Nikitin et al., 2012). They are often characterized by the serial correlation coefficient (SCC)

$\begin{array}{cc}{\rho }_k=\frac{\langle \left(T_i-\langle T_i\rangle \right)\left(T_{i+k}-\langle T_{i+k}\rangle \right)\rangle }{\langle {(T_i-\langle T_i\rangle )}^2\rangle },& (1)\end{array}$

where T_i and T_{i + k} are two ISIs lagged by an integer k and 〈·〉 denotes ensemble averaging. ISI correlations can be induced via correlated input to the neural dynamics, e.g. in the form of external colored noise (Middleton et al., 2003; Lindner, 2004), intrinsic noise from ion channels with slow kinetics (Fisch et al., 2012), or stochastic narrow-band input (Neiman and Russell, 2001, 2005; Bauermeister et al., 2013).

Another ubiquitous mechanism for ISI correlations are slow feedback processes mediating spike-frequency adaptation (Chacron et al., 2000; Liu and Wang, 2001; Benda et al., 2005)—a phenomenon describing the reduced neuronal response to slowly changing stimuli (Benda and Herz, 2003; Gabbiani and Krapp, 2006). In the stationary state, these adaptation mechanisms are typically associated with short-range correlations with a negative SCC at lag k = 1 and a reduced Fano factor as demonstrated by several numerical (Geisler and Goldberg, 1966; Wang, 1998; Liu and Wang, 2001; Benda et al., 2010) and analytical studies (Schwalger et al., 2010; Schwalger and Lindner, 2010; Farkhooi et al., 2011; Urdapilleta, 2011). The correlation structure of adapting neurons can show qualitatively different patterns, ranging from monotonically decaying correlations to damped oscillations when plotted as a function of the lag (Ratnam and Nelson, 2000). Because ISI correlations shape spectral measures (Chacron et al., 2004), they bear implications for neural computation in general. However, a simple theory that predicts and explains possible correlation patterns is still lacking.

In this article, we present a relation between the ISI correlation coefficient ρ_k and a basic characteristics of nonlinear neural dynamics, the phase-response curve (PRC). The PRC quantifies the advance (or delay) of the next spike caused by a small depolarizing current applied at a certain time after the last spike (Ermentrout, 1996). For neurons which integrate up their input (integrator neurons), the PRC is positive at all times (type I PRC) whereas neurons, which show subthreshold resonances (resonator neurons), possess a PRC that is partly negative (type II PRC) (Ermentrout, 1996; Izhikevich, 2005; Ermentrout and Terman, 2010). Below we show that resonator neurons possess a richer repertoire of correlation patterns than integrator neurons do.

2. Results

2.1. Model

Spike frequency adaptation can be modeled by Hodgkin–Huxley type neurons with a depolarization-activated adaptation current (Wang, 1998; Ermentrout et al., 2001; Benda and Herz, 2003). However, the spiking of such conductance-based models can in many instances be approximated by simpler multi-dimensional integrate-fire (IF) models that are equipped with a spike-triggered adaptation current (Treves, 1993; Izhikevich, 2003; Brette and Gerstner, 2005); adapting IF models perform excellently in predicting spike times of real cells under noisy stimulation (Gerstner and Naud, 2009). Here, we consider a stochastic nonlinear multi-dimensional IF model for the membrane potential v, N auxiliary variables w_j (j = 1, …, N) and a spike-triggered adaption current a(t):

$\begin{array}{cc}\dot{v}=f_0\text{}\left(v,w\right)+\mu -a+\xi (t),& (2\text{a})\end{array}$

$\begin{array}{cc}{\dot{w}}_j=f_j(v,w),& (2\text{b)}\end{array}$

$\begin{array}{cc}{\tau }_a\dot{a}=-a+{\tau }_a\Delta {\displaystyle \sum _i\delta }\left(t-t_i\right).& (2\text{c)}\end{array}$

The membrane potential v(t) is subject to weak Gaussian noise ξ(t) with 〈ξ(t)ξ(t′)〉 = 2Dδ(t − t′) and noise intensity D. The dynamics is complemented by a spike-and-reset mechanism: whenever v_(T) reaches a threshold υ(t), a spike is registered at time t_i = t and v(t) and w(t) = [w₁(t), …, w_N(t)]^T are reset to v(t⁺_i) = 0 and w(t⁺_i) = w_r (where t⁺_i denotes the right-sided limit t → t_i + 0). At the same time, a(t) suffers a jump by Δ ≥ 0 as seen from Equation (2c), which resembles high-threshold adaptation currents (Wang, 1998; Liu and Wang, 2001). The constant input current μ is assumed to be sufficiently large to ensure ongoing spiking even in the absence of noise. Note that the model is non-dimensionalized by measuring time in units of the membrane time constant τm~ 10 ms and voltage in units of the distance between reset and spike-initiating potential (a typical value is 15 mV). In particular, the adaptation time constant τ_a is measured relative to τm and the unit of the firing rate is τ⁻¹_m ~ 100 Hz.

An important special case, the adaptive exponential integrate-and-fire model (Brette and Gerstner, 2005) with purely spike-triggered adaptation and a white noise current with constant mean is illustrated in Figure 1. It assumes an exponential nonlinearity f₀(v) = −γv + γΔ_T exp[(v − 1)/Δ_T] (Fourcaud-Trocmé et al., 2003; Badel et al., 2008) and corresponds to N = 0. Time courses of v(t) and a(t) are shown in Figures 1A1,B1 for two distinct correlation patterns possible in this model. The ISIs T_i = t_i − t_{i − 1} are obtained as differences between subsequent spiking times t_i. The sequence T_i, T_{i + 1}, T_{i + 2} displays patterns of short-long-long (Figure 1A1) and short-long-short (Figure 1B1), corresponding to a negative SCC, which decays monotonically with the lag k (Figure 1A3) or to an SCC oscillating with k (Figure 1B3). In the following, we develop a theory to analyze these and other correlation patterns possible in multi-dimensional adapting IF models.

FIGURE 1

Figure 1. Correlation patterns in the adaptive exponential IF model with τ_a = 10, γ = 1, Δ_T = 0.1, υ_T = 2, D = 0.1. Adaptation is weak (Δ = 1, μ = 15) in (A) and strong (Δ = 10, μ = 80) in (B). Membrane voltage v(t) and adaptation variable a(t) with ISI sequences {T_i} and peak adaptation values {a_i} are shown in (A1,B1); time is in units of the membrane time constant τm. Colored pieces of trajectories in the phase plane (v,a) in (A2,B2) correspond to the respective colors in (A1,B1). The deterministic limit cycle (LC), determined by the initial (post-spike) values v = 0, a = a*, is indicated by a thick black line. For weak adaptation (A2) a short ISI T_i causes positive deviations δa_i = a_i − a* and δa_{i + 1} = a_{i + 1} − a* of peak values leading to long ISIs T_{i + 1} and T_{i + 2} and, hence, to a negative ISI correlation at all lags (A3). Because of the qualitatively different limit cycle for strong adaptation (B2), deviations δa_i and δa_{i + 1} differ in sign, yielding an oscillatory correlation pattern (B3).

2.2. General Theory

In our model Equation (2), a(t) is the only variable that keeps a memory of the previous spike times thereby inducing correlations between ISIs. Over one ISI the time course of adaptation is an exponential decay, relating two adjacent peak values a_i = a(t⁺_i) and a_{i + 1} = a(t⁺_{i + 1}) by

$\begin{array}{cc}{a}_{i+1}=a_i{e}^{-T_{i+1}/{\tau }_a}+\Delta & (3)\end{array}$

(Figures A1,B1). We assume that in the deterministic case (D = 0) our model has a finite period T* (i.e., the model operates in the tonically firing regime) and, hence, for D = 0 the map (3) has a stable fixed point

$\begin{array}{cc}{a}^{\ast }=\Delta /\left[1-\exp \left(-T^{\ast }/{\tau }_a\right)\right].& (4)\end{array}$

The asymptotic deterministic dynamics can be interpreted as a limit-cycle like motion in the phase space from the reset point to the threshold and back by the instantaneous reset [cf. Figures 1A2,B2].

Weak noise will cause small deviations in the period δT_i = T_i − T*≈ T_i − 〈T_i〉 that are mutually correlated with coefficient ρ_k = 〈δT_iδT_{i + k}〉/〈δT²_i〉. The peak adaptation values, however, also fluctuate, δa_i = a_i − a*, and both deviations are related by linearizing Equation (3):

$\begin{array}{cc}\delta T_{i+1}=\frac{{\tau }_a}{a^{\ast }}\left(\delta a_i-e^{T^{\ast }/{\tau }_a}\delta a_{i+1}\right).& (5)\end{array}$

A second relation between the small deviations can be gained by considering how a small perturbation in the voltage dynamics affects the length of the period. This effect is captured by the infinitesimal phase response curve (PRC), Z(t), t ∈ (0, T*) (Izhikevich, 2005; Ermentrout and Terman, 2010) (see Section 4 for the precise definition). During the interval T_{i + 1}, the voltage dynamics in Equation (2a) can be written as $\dot{v}=f_0(v,w)+\mu -(a^{\ast }+\delta a_i)e^{-(t-t_i)/{\tau }_a}+\xi (t)$. Compared to the deterministic limit cycle, the dynamics is perturbed by the weak noise and the small deviation in the adaptation δa_ie^{−(t − t_i)/τ_a} yielding in linear response

$\begin{array}{cc}\delta T_{i+1}={\displaystyle {\int }_0^{T^{\ast }}\text{d}t}Z(t)\left(\delta a_i{e}^{-\frac{t}{{\tau }_a}}-\xi (t_i+t)\right).& (6)\end{array}$

Combining Equations (5), (6) we obtain the stochastic map

$\begin{array}{cc}\delta a_{i+1}=\alpha \vartheta \delta a_i+{\Xi }_i,& (7)\end{array}$

where ${\Xi }_i=\frac{\alpha a^{*}}{{\tau }_a}{\displaystyle {\int }_0^{\infty }\text{d}t}Z(t)\xi (t_i+t)$ are uncorrelated Gaussian random numbers and

$\begin{array}{cc}\alpha =e^{-T^{\ast }/{\tau }_a},\vartheta =1-\frac{a^{\ast }}{{\tau }_a}{\displaystyle {\int }_0^{T^{\ast }}\text{d}t}Z(t)e^{-\frac{t}{{\tau }_a}}.& (8)\end{array}$

Note that local stability of the fixed point a* requires that |αϑ| < 1. The covariance c_k = 〈δa_iδa_{i + k}〉 of the auto-regressive process Equation (7) can be calculated by elementary means and using Equation (5) we obtain for k ≥ 1:

$\begin{array}{cc}{\rho }_k=-A(1-\vartheta ){\left(\alpha \vartheta \right)}^{k-1},A=\frac{\alpha (1-{\alpha }^2\vartheta )}{1+{\alpha }^2-2{\alpha }^2\vartheta }.& (9)\end{array}$

In order to compute α and ϑ via Equation (8), we have to calculate T* and Z(t) (a* then follows from Equation (4)), which can be done for simple systems analytically.

Our main result, Equations (8), (9), allows to draw a number of general conclusions. It shows that the SCC is always a geometric sequence with respect to the lag k that can generate qualitatively different correlation patterns depending on the value of ϑ and thus on PRC and adaptation current. Because |αϑ| < 1 and 0 < α < 1, the prefactor A in Equation (9) is always positive. Consequently, ρ₁ is negative for ϑ < 1 and positive for ϑ > 1. Looking at Equation (8), we find that a positive PRC inevitably yields ϑ < 1. This implies that adapting neurons with type I PRC possess negative correlations between adjacent ISIs. Intuitively, a short ISI causes in the following on average a higher inhibitory adaptation during the subsequent ISI. Such an inhibitory current always enlarges the ISI in type I neurons—hence, a short ISI is followed by a long ISI.

The sign of higher lags is determined by the base of the power: for ϑ > 0 correlations decay monotonically, whereas for ϑ < 0 the SCC oscillates. Two special cases are ϑ = 0 with a negative correlation at lag 1 and vanishing correlations at all higher lags and ϑ = 1 where all correlations vanish. Overall, we find five basic patterns corresponding to the cases −α⁻¹ < ϑ < 0, ϑ = 0, 0 < ϑ < 1, ϑ = 1 and 1 < ϑ < α⁻¹. These basic patterns cover all interval correlations discussed in previous theoretical studies (Schwalger and Lindner, 2010; Urdapilleta, 2011). Our geometric formula generalizes the theory for the perfect IF model with adaptation (Schwalger et al., 2010) to more realistic, nonlinear multi-dimensional IF models with adaptation.

The cumulative effect of the correlations can be described by the sum over all ρ_k, which determines the long-time limit of the Fano factor and the low-frequency limit of the spike train power spectrum (for a definition of these quantities, see Section 4.2). Evaluating the geometric series yields

$\begin{array}{cc}{\displaystyle \sum _{k=1}^{\infty }{\rho }_k}=-\frac{A\left(1-\vartheta \right)}{1-\alpha \vartheta }.& (10)\end{array}$

This shows that adaptation in neurons with type I resetting (ϑ < 1) leads to a negative summed correlation and hence a reduced long-term variability. Furthermore, at high firing rates achieved by a strong input current μ, the sum in Equation (10) can be approximated by

$\begin{array}{cc}{\displaystyle \sum _{k=1}^{\infty }{\rho }_k}\simeq -\frac{1}{2}+\frac{1/2}{{\left(1+\Delta {\tau }_a/v_{\text{T}}\right)}^2},T^{\ast }\ll {\tau }_a.& (11)\end{array}$

In particular, for strong adaptation (Δτ_a ≫ v_T) the sum is only slightly larger than −1/2. Note that by virtue of the fundamental relation ${\lim }_{t\to \infty }F(t)=C_{\text{V}}^2\left(1+2{\displaystyle {\sum }_{k=1}^{\infty }{\rho }_k}\right)$ (Cox and Lewis, 1966) (see Section 4.2), the smallest possible value for the sum is −1/2 in order to ensure the non-negativity of the Fano factor F(t). At this minimal value the long-term variability as expressed by the Fano factor vanishes even for a non-vanishing ISI variability as quantified by the coefficient of variation C_V. The latter quantity can also be estimated using the weak-noise theory: From Equation (7) one can calculate the variance of a_i and using Equation (5) an approximation for C²_V ≈ 〈δT²_i〉/T^*2 can be obtained as follows:

$\begin{array}{cc}{C}_{\text{V}}^2=2D\frac{1+{\alpha }^2-2{\alpha }^2\vartheta }{[1-{\left(\alpha \vartheta \right)}^2]T^{\ast }{}^2}{\displaystyle {\int }_0^{T^{\ast }}\text{d}t{[Z(t)]}^2}.& (12)\end{array}$

2.3. One-Dimensional IF Models with Adaptation

In the simplest case (N = 0, f₀(v, w) = f(v)) the PRC reads

$\begin{array}{cc}Z(t)=Z(T^{\ast })\exp [{\displaystyle {\int }_t^{T^{\ast }}\text{d}{t}^{\prime }}{f}^{\prime }(v_0(t^{\prime }))],& (13)\end{array}$

where v₀(t) is the limit cycle solution and Z(T*) = [f(v_T) + μ − a* + Δ]⁻¹ is the inverse of the velocity ${\dot{v}}_0(T^{\ast })$ at the threshold, which is always positive. Thus, the PRC is positive for all t ∈ (0, T*), i.e., one-dimensional IF models show type I behavior. From our general considerations, this implies a negative SCC at lag k = 1. The sign of the correlations at higher lags can be inferred from the sign of ϑ, for which one can show (Section 4) that

$\begin{array}{cc}\vartheta =\left(f(0)+\mu -a^{\ast }\right)Z(0).& (14)\end{array}$

Because Z(0) > 0, the sign of ϑ is determined by the sign of f(0) + μ − a*. For weak adaptation such that a* < f(0) + μ (achieved by a sufficiently small value of Δ or τ_a, Figure 1A), we will have ϑ > 0 and a negative correlation at all lags (Figure 1A3). In this case, a short ISI occurring by fluctuation will cause a positive deviation δa_i (Figure 1A2, green arrow). Geometrically, it is plausible that such a positive deviation causes a likewise positive deviation δa_{i + 1} in the subsequent cycle (Figure 1A2, red arrow). Because a positive deviation is associated with a long ISI, the initial short ISI is on average followed by longer ISIs.

In marked contrast, for strong adaptation such that a* > f(0) + μ (achieved by a sufficiently large value of Δ or τ_a), ϑ becomes negative and hence the SCC's sign alternates with the lag. This alternation of the sign can be understood by means of the phase plane. Let us again consider a positive deviation δa_i due to a short preceding ISI (Figure 1B2, green arrow). Because ${\dot{v}}_0(0)=f(0)+\mu -a^{\ast }<0$, the neuron is reset above the v-nullcline and hence hyperpolarizes at the beginning of the interval, i.e., the trajectory makes a detour into the region of negative voltage (corresponding to a “broad reset” in Naud et al. (2008)). A positive deviation δa_i leads to a larger detour (green trajectory) causing a sign inversion and hence a negative deviation δa_{i + 1} (Figure 1B2, red arrow). Because a positive (negative) deviation corresponds on average to a long (short) ISI, the alternation of δa_i also entails an alternation of the ISI correlations. Thus, the distinction between monotonic and alternating patterns relates to a qualitative distinction of the voltage trace after resetting [cf. “sharp” vs. “broad” resets in Naud et al. (2008)].

As demonstrated in Figures 1A3,B3, our theory works well for the adapting exponential integrate-and-fire model. We next demonstrate the validity of our approach over a broad range of firing rates (Figure 2) for another important 1D model, the adapting leaky integrate-and-fire model (Treves, 1993) for which f(v) = −γv and

$\begin{array}{cc}Z(t)=\exp \left[\gamma (t-T^{\ast })\right]/\left(\mu -\gamma v_{\text{T}}-a^{\ast }+\Delta \right)& (15)\end{array}$

(here T* has still to be determined from a transcendental equation). Changing the firing rate by varying the input current μ, we find a good agreement for the first two correlation coefficients and the sum of all ρ_k; the approximation of the CV shows deviations from simulation results when the input current μ becomes small (approaching the fluctuation-driven regime). In accordance with previous findings (Wang, 1998; Liu and Wang, 2001; Benda et al., 2010; Nesse et al., 2010; Schwalger et al., 2010; Schwalger and Lindner, 2010; Urdapilleta, 2011), the first correlation coefficient ρ₁ displays a minimum corresponding to strong anti-correlations between adjacent intervals. The correlations at lag 2 can be positive for a finite range of firing rates if the adaptation strength is sufficiently large (Figure 2B), whereas for moderate adaptation we find a negative ρ₂ at all firing rates (Figure 2A). In both cases, however, the sum of SCCs approaches a value close to −1/2 for high firing rates as predicted by Equation (11) (Figure 2, bottom). This is strikingly similar to experimental data from weakly electric fish, in which some electro-receptors display a monotonically decaying SCC and some show an oscillatory SCC (Ratnam and Nelson, 2000) but all cells exhibit a sum close to −1/2 (Ratnam and Goense, 2004). Finally, Figure 2 reveals a local maximum of the CV for some suprathreshold current μ—an effect that has been described by Nesse et al. (2008).

FIGURE 2

Figure 2. ISI correlations and coefficient of variation (CV) of the adapting LIF model vs. firing rate 1/〈T_i〉≈ 1/T*, where the rate is varied by increasing μ. The gray-shaded area corresponds to the fluctuation-driven regime (μ < γv_T), where the assumptions of the theory do not hold. The panels display (from top to bottom) ρ₁, ρ₂, the sum ${\displaystyle {\sum }_{k=1}^m{\rho }_k}$ and the CV for simulation (circles, m = 100) and theory (solid lines, m → ∞). (A) Moderate adaptation: Δ = 1, (B) strong adaptation: Δ = 10. Both: γ = 1, τ_a = 10, D = 0.1, v_T = 1. Note that the firing rate is given in units of the inverse membrane time constant τ⁻¹_m.

2.4. Generalized Integrate-and-Fire Model with Adaptation

Different correlation patterns become possible if we consider a type II PRC, which is by definition partly negative and can lead to a negative value of the integral in Equation (8), and hence to ϑ ≥ 1. This corresponds to a non-negative SCC at lag 1, which is infeasible in the one-dimensional case. To test the prediction ρ₁ ≥ 0, we study the generalized integrate-and-fire (GIF) model (Brunel et al., 2003) with spike-triggered adaptation. This model is defined by f₀(v, w) = −γv − βw and f₁(v, w) = (v − w)/τ_w. Using the method described in Section 4, the PRC is obtained as

$\begin{array}{cc}Z(t)=\frac{e^{\frac{\nu }{2}(t-T^{\ast })}\left[\cos (\Omega (t-T^{\ast }))-\frac{1-{\tau }_w\gamma }{2{\tau }_w\Omega }\sin \left(\Omega (t-T^{\ast })\right)\right]}{\mu -\gamma v_{\text{T}}-\beta w_0(T^{\ast })-a^{\ast }+\Delta }& (16)\end{array}$

where $\nu =\gamma +1/{\tau }_w,\Omega =\sqrt{\frac{\beta +\gamma }{{\tau }_w}-\frac{{\nu }^2}{4}}$ and w₀(t) is one component of the deterministic limit-cycle solution [v₀(t), w₀(t), a₀(t)] that we calculated numerically.

In Figure 3B we demonstrate that all possible correlation patterns can be realized in the GIF model and that the predicted SCCs agree quantitatively well in theory and model simulations (for comparison, see the SCC for the LIF in Figure 3A). To each distinct pattern belongs a range of ϑ (Figure 3, left), determined by the area under the weighted PRC $\tilde{Z}(t)=\frac{a^{\ast }}{{\tau }_a}{e}^{-\frac{t}{{\tau }_a}}Z(t)$. The function $\tilde{Z}(t)$ (left column in Figures 3A,B) illustrates, why an adapting GIF neuron can show vanishing (Figure 3Biv) or even purely positive ISI correlations (Figure 3Bv). In case of type II resetting, inhibitory input can shorten the ISI because of the negative part in the PRC; here inhibition acts like an excitatory input. Consequently, a short ISI will induce a stronger inhibition (adaptation) that now causes a likewise short interval and results thus in a positive correlation between adjacent ISIs. Also, the shortening effect of the adaption current in the early negative phase of the PRC can be exactly balanced by the delaying effect of the late positive phase of the PRC (pseudo-renewal case, in which the area under $\tilde{Z}$ is zero).

FIGURE 3

Figure 3. Possible correlation patterns and corresponding PRCs (solid lines: theory, symbols: simulations of Equation (2)). For the adapting LIF model (A), ϑ < 1 and only three qualitative different cases are possible. The adapting GIF model (B) exhibits the full repertoire of correlation patterns because the PRC can be partly negative and ϑ can attain values from its entire physically meaningful interval [−1/α, 1/α]. The value of ϑ and hence the type of correlation pattern is set by the integral over the weighted PRC $\tilde{Z}(t)=Z(t)e^{-\frac{t}{{\tau }_a}}\frac{a^{\ast }{T}^{\ast }}{{\tau }_a}$, shown in left panels. LIF parameters: D = 0.1, τ_a = 2, (i) μ = 20, Δ = 10, (ii) μ = 20, Δ = 4.47, (iii) μ = 5, Δ = 1. GIF parameters: (i) μ = 10, β = 3, τ_a = 10. (ii) μ = 11.75, β = 3, τ_a = 10. (iii) μ = 20, β = 1.5, τ_a = 10. (iv) μ = 2.12, β = 1.5, τ_a = 1, Δ = 10. (v) μ = 1.5, β = 1.5, τ_a = 1, Δ = 9 D = 10⁻⁵. Unless stated otherwise, γ = 1, Δ = 1, τ_w = 1.5, D = 10⁻⁴, w_r = 0.

3. Discussion

We have found a general relation between two experimentally accessible characteristics: the serial interval correlations and the phase response curve of a noisy neuron with spike-triggered adaptation. The theory predicts distinct correlation patterns like short-range negative and oscillatory correlations that have been observed in experiments (Ratnam and Nelson, 2000; Nawrot et al., 2007) and in simulation studies of adapting neurons (Chacron et al., 2000; Liu and Wang, 2001).

Beyond negative and oscillatory correlations, we have found, however, that resonator neurons with spike-frequency adaptation can exhibit purely positive ISI correlations or a pseudo-renewal process with uncorrelated intervals. Adaptation currents that are commonly associated with negative ISI correlations (Wang, 1998; Chacron et al., 2001; Liu and Wang, 2001; Chacron et al., 2003; Benda et al., 2010; Nesse et al., 2010) can thus induce a rich repertoire of correlation patterns. Despite the multitude of patterns, there is a universal limit for the cumulative correlations at high firing rates [cf. Equation (11)], which shows that the long-term variability of the spike train is in this limit always reduced in agreement with experimental studies (Ratnam and Goense, 2004).

Our analytical results apply to arbitrary adaptation strength and time scale but require that (1) the noise is weak and white, (2) the deterministic dynamics shows periodic firing with equal ISIs (i.e., a limit-cycle exists) and (3) the adaptation current is purely spike-triggered with (4) a single exponential decay time. Regarding the weak-noise assumption, we found from numerical simulations quantitative agreement with our theory for values of the coefficient of variation (CV) up to 0.4, which is, for instance, typical for neurons in the sensory periphery (Ratnam and Nelson, 2000; Neiman and Russell, 2004; Vogel et al., 2005). This holds even in the subthreshold regime at low CVs, where the deterministic system does not follow a limit cycle. In this case, T* has to be replaced by the mean ISI. Moreover, we found qualitative agreement even for moderately strong noise with values of the CV up to 0.8, which is typical for cortical non-bursting neurons in vivo (e.g. Figure 3 in Softky and Koch (1993)).

In the absence of a deterministic limit-cycle, i.e., in the fluctuation-driven regime at high CVs, different mathematical approaches have to be employed, such as those based on a hazard-function formalism (Muller et al., 2007; Nesse et al., 2010; Schwalger and Lindner, 2010; Farkhooi et al., 2011). Furthermore, for some parameter sets, we also observed repeat periods of the deterministic system that involved multiple ISIs corresponding to a periodic ISI sequence with T_i = T_{i + n}, where the smallest period is n ≥ 2. Such cases can realize bursting (Naud et al., 2008), which we did not consider in the present study. However, we expect that these parameter regimes yield interesting correlation patterns because already in the noiseless case a periodic ISI sequence exhibits correlations between ISIs.

Regarding the last two assumptions, it seems that the analytical derivation cannot be easily extended to the cases of adaptation currents activated by the subthreshold membrane potential (“subthreshold adaptation” Ermentrout et al., 2001; Brette and Gerstner, 2005; Prescott and Sejnowski, 2008; Deemyad et al., 2012) and multiple-time-scale adaptation (Pozzorini et al., 2013). Ermentrout et al. (2001) have shown that the inclusion of subthreshold adaptation can lead to type II PRCs, which according to our theory could qualitatively change the correlation patterns. An adaptation dynamics depending on the subthreshold membrane potential also involves a fluctuating component because v is noisy. According to Schwalger et al. (2010), this stochasticity could contribute positive correlations. The combined effect of spike-triggered, subthreshold and stochastic adaptation currents on the sign of the SCC is not clear.

The important cases of the fluctuation-driven regime and multiple-time-scale adaptation have been recently analyzed with respect to the first-order spiking statistics including the stationary firing rate as well as the mean response to time-dependent stimuli (Richardson, 2009; Naud and Gerstner, 2012). The second-order statistics, which describes the fluctuations of the spike train (“neural variability,” cf. Section 4.2) and which limits the information transmission capabilities of neurons, is however still poorly understood theoretically in these cases. How adaptation shapes second-order statistics in the cases of multiple adaptation time scales, fluctuation-driven spiking and sub-threshold adaptation is an interesting topic for future investigations.

As an outlook we sketch, how our theory could be used to constrain unknown physiological parameters by measured SCCs and PRCs. For instance, from the mean ISI we can estimate T* = 〈T〉. Furthermore, knowing ρ₁ = −A(α, ϑ)(1 − ϑ) as well as the ratio ρ₂/ρ₁ = αϑ one can eliminate ϑ and solve for α. This allows to estimate the unknown adaptation time constant τ_a = −T*/ln α and the amplitude of the adaptation current

$\begin{array}{cc}{a}^{\ast }=\frac{{\tau }_a}{\alpha }\left(\alpha -\frac{{\rho }_2}{{\rho }_1}\right)/{\displaystyle {\int }_0^{T^{\ast }}\text{d}t}Z(t)e^{-\frac{t}{{\tau }_a}}.& (17)\end{array}$

Although experimental PRCs are notoriously noisy (Izhikevich, 2005), the integral over Z(t) determining our estimate of a* is less error-prone. Combining our approach with advanced estimation methods for the PRC (Galán et al., 2005), may thus provide an alternative access to hidden physiological parameters using only spike time statistics.

4. Materials and Methods

4.1. Phase-Response Curves of Adapting IF Models

We use the phase-response curve Z(t′) to characterize the shift of the next spike following a small current pulse applied at a given “phase” t' ∈ [0, T*] of an ISI. More precisely, let us assume that the last spike occurred at time t₀ = 0. Then, the next spike time t₁ of the perturbed limit cycle dynamics $\dot{v}=f_0(v,w)+\mu -a+ϵ\delta (t-t^{\prime })$, v(0) = 0, w(0) = w_r, a(0) = a*, 0 < t′ ≤ T* will be shifted by some amount δT(t′, ϵ) = t₁ − T*. The infinitesimal PRC can be defined as the limit

$\begin{array}{cc}Z(t^{\prime })=-\underset{ϵ\to 0}{\lim }\frac{\delta T(t^{\prime },ϵ)}{ϵ},& (18)\end{array}$

where the sign has been chosen such that a spike advance (δT < 0) due to a positive stimulation (ϵ > 0) leads to a positive PRC. The definition of Z(t) by the shift of the next spike differs from the PRC that describes the asymptotic spike shift but is equivalent to the so-called “first-order PRC,” which is often measured in experiments (Netoff et al., 2012).

4.1.1. Adjoint equation and boundary conditions

The PRC can be computed using the adjoint method (see e.g. Ermentrout and Terman (2010)). To this end, the dynamics is linearized about the T*-periodic limit cycle solution y₀(t) = [v₀(t), w₀(t), a₀(t)]. The linearized limit-cycle dynamics y(t) = y₀(t) + δy(t) corresponding to Equation (2) is given by

$\begin{array}{cc}\stackrel{.}{\delta y}=A(t)\delta y& (19)\end{array}$

with the Jacobian matrix

$\begin{array}{cc}A(t)=\left(\begin{array}{l}\frac{\partial f_0}{\partial v}\frac{\partial f_0}{\partial w_1}\dots \frac{\partial f_0}{\partial w_N}-1\\ {\tau }_1^{-1}\frac{\partial f_1}{\partial v}{\tau }_1^{-1}\frac{\partial f_1}{\partial w_1}\dots {\tau }_1^{-1}\frac{\partial f_1}{\partial w_N}0\\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \\ {\tau }_N^{-1}\frac{\partial f_N}{\partial v}{\tau }_N^{-1}\frac{\partial f_N}{\partial w_1}\dots {\tau }_N^{-1}\frac{\partial f_N}{\partial w_N}0\\ 0\cdot \cdot \cdot 0-{\tau }_a^{-1}\end{array}\right)& (20)\end{array}$

evaluated at v = v₀(t), w = w₀(t). The linear response of the ISI to perturbations of the limit-cycle dynamics in an arbitrary direction is given by the vector Z(t) = [Z(t), Z_w1(t), …, Z_wN(t), Z_a(t)]^T, where the first component is equal to the PRC defined above. This vector satisfies the adjoint equation Ż = −A^T Z (A^T denotes the transpose of A) with the normalization condition ${\dot{v}}_0(t)Z(t)+{\dot{w}}_0(t)Z_w(t)+{\dot{a}}_0(t)Z_a(t)=1$. The remaining N + 1 boundary conditions are obtained by the following consideration: On the limit cycle Γ, a phase ϕ: Γ → [0, T*] can be introduced in the usual way by inverting the map t ↦ y₀(t) and setting ϕ = t. Because we are interested in the shift of the next spike, it is useful to define the isochrons (sets of equal phase) as the sets of all points in phase space that will lead to the same first spike time. Put differently, phase points belonging to the same isochron will have their first threshold crossing in synchrony. As a consequence, the threshold hyperplane defined by the condition v = v_T is a special isochron corresponding to the phase ϕ = T*. Note that this definition of the phase implies that the reset line defined by the condition v = 0, w = w_r does generally not correspond to ϕ = 0 but to positive phases if a < a* and negative phases if a > a*. Thus, off-limit-cycle trajectories suffer a phase jump upon reset. Close to the threshold, the isochrons are parallel to the threshold, and thus, a perturbation perpendicular to the v-direction does not change the phase. This insensitivity implies the boundary conditions Z_w1(T*) = … = Z_wN(T*) = Z_a(T*) = 0. Note that a definition of the PRC based on the asymptotic spike shift would require periodic boundary conditions (Ladenbauer et al., 2012).

From the above considerations, it becomes clear that the PRC Z(t) can be computed for t ∈ [0, T*] by solving the system

$\begin{array}{cc}\left(\begin{array}{c}\dot{Z}\\ {\dot{Z}}_{w_1}\\ ⋮\\ {\dot{Z}}_{w_N}\end{array}\right)=-\left(\begin{array}{cccc}\frac{\partial f_0}{\partial v}& {\tau }_1^{-1}\frac{\partial f_1}{\partial v}& \dots & {\tau }_N^{-1}\frac{\partial f_N}{\partial v}\\ \frac{\partial f_0}{\partial w_1}& {\tau }_1^{-1}\frac{\partial f_1}{\partial w_1}& \dots & {\tau }_N^{-1}\frac{\partial f_N}{\partial w_1}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f_0}{\partial w_N}& {\tau }_1^{-1}\frac{\partial f_1}{\partial w_1}& \dots & {\tau }_N^{-1}\frac{\partial f_N}{\partial w_N}\end{array}\right)\left(\begin{array}{c}Z\\ Z_{w_1}\\ ⋮\\ Z_{w_N}\end{array}\right)& (21)\end{array}$

subject to the boundary conditions

$\begin{array}{cc}Z(T^{\ast })=\frac{1}{{\dot{v}}_0(T^{\ast })}=\frac{1}{f_0(v_{\text{T}},w_0(T^{\ast }))+\mu -a^{\ast }+\Delta },& (22)\end{array}$

$\begin{array}{cc}{Z}_{w_k}(T^{\ast })=0,k=1,\dots ,N.& (23)\end{array}$

The PRC with respect to a is determined by

$\begin{array}{cc}{\dot{Z}}_a=\frac{1}{{\tau }_a}{Z}_a+Z(t),Z_a(T^{\ast })=0.& (24)\end{array}$

The matrix in Equation (21) is again evaluated on the limit cycle at v = v₀(t), w = w₀(t) and is therefore time-dependent. An analytic solution of Equation (21) is possible for one-dimensional models with adaptation (N = 0) or general linear IF models although in most cases the deterministic period T* still has to be computed numerically.

4.1.2. One-dimensional case

In the case N = 0, the PRC satisfies the equation Ż = −f′(v₀)Z with boundary condition (22). The solution is given by Equation (13). In order to prove Equation (14), we compute Z_a(t) from Equation (24) yielding

$Z_a(t)=e^{\frac{t}{{\tau }_a}}\left(Z_a(0)+{\displaystyle {\int }_0^{t}Z(t^{\prime })e^{-\frac{t^{\prime }}{{\tau }_a}}}\text{ d}{t}^{\prime }\right).$

Evaluation of this expression for t = T* leads to $\vartheta =1+\frac{a^{*}}{{\tau }_a}{Z}_a(0)$. Finally, using the normalization condition $(f(0)+\mu -a^{\ast })Z(0)-\frac{a^{\ast }}{{\tau }_a}{Z}_a(0)=1$ yields Equation (14).

4.2. Relation Between Second-Order Statistics of Spike Count, Spike Train and Interspike Intervals

A stationary sequence of spike times {…, t_{i − 1}, t_i, t_{i + 1}, …} is often characterized by the statistics of the spike train x(t) = ∑_iδ(t − t_i), the spike count $N(t)={\displaystyle {\int }_0^t\text{d}{t}^{\prime }}x(t^{\prime })$ or the sequence of ISIs {T_i = t_i − t_{i − 1}}. In particular, neural variability can be quantified by the second-order statistics of these different descriptions as, for instance, the spike train power spectrum

$\begin{array}{cc}S(f)={\displaystyle \int \text{d}}\tau e^{2\pi if\tau }\langle x(t)x(t+\tau )\rangle ,& (25)\end{array}$

the Fano factor

$\begin{array}{cc}F(t)=\frac{\langle N{(t)}^2\rangle -{\langle N(t)\rangle }^2}{\langle N(t)\rangle },& (26)\end{array}$

and the coefficient of variation $C_{\text{V}}=\sqrt{\langle {(T_i-\langle T_i\rangle )}^2\rangle }/\langle T_i\rangle$ and SCC ρ_k as defined in Equation (1). These statistics are connected by the fundamental relationship (Cox and Lewis, 1966) (see also (van Vreeswijk, 2010))

$\begin{array}{cc}\underset{t\to \infty }{\lim }F(t)=\langle T_i\rangle \underset{f\to 0}{\lim }S(f)=C_{\text{V}}^2\left(1+2{\displaystyle \sum _{k=1}^{\infty }{\rho }_k}\right).& (27)\end{array}$

It shows that the summed SCC has a strong impact on the long-term variability of the spike train. In particular, a negative sum yields a more regular spike train on long time scales than a renewal process with the same C_V.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by Bundesministerium für Bildung und Forschung grant 01GQ1001A.

References

Avila-Akerberg, O., and Chacron, M. J. (2011). Nonrenewal spike train statistics: causes and functional consequences on neural coding. Exp. Brain Res. 210, 353–371. doi: 10.1007/s00221-011-2553-y

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Badel, L. Lefort, S. Brette, R. Petersen, C. C., Gerstner, W., and Richardson., M. J. (2008). Dynamic i-v curves are reliable predictors of naturalistic pyramidal-neuron voltage traces. J. Neurophysiol. 99, 656–666. doi: 10.1152/jn.01107.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bauermeister, C. Schwalger, T. Russell, D. Neiman, A., and Lindner., B. (2013). Characteristic effects of stochastic oscillatory forcing on neural firing statistics: theory and application to paddlefish electroreceptor afferents. PLoS Comput. Biol. 9:e1003170. doi: 10.1371/journal.pcbi.1003170

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Benda, J., and Herz, A. V. M. (2003). A universal model for spike-frequency adaptation. Neural Comp. 15, 2523. doi: 10.1162/089976603322385063

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Benda, J., Longtin, A., and Maler, L. (2005). Spike-Frequency adaptation separates transient communication signals from background oscillations. J. Neurosci. 25, 2312–2321. doi: 10.1523/JNEUROSCI.4795-04.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Benda, J., Maler, L., and Longtin, A. (2010). Linear versus nonlinear signal transmission in neuron models with adaptation currents or dynamic thresholds. J. Neurophysiol. 104, 2806–2820. doi: 10.1152/jn.00240.2010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Brette, R., and Gerstner, W. (2005). Adaptive exponential integrate-and-fire model as an effective description of neuronal activity. J. Neurophysiol. 94, 3637–3642. doi: 10.1152/jn.00686.2005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Brunel, N., Hakim, V., and Richardson, M. J. (2003). Firing-rate resonance in a generalized integrate-and-fire neuron with subthreshold resonance. Phys. Rev. E 67(5 Pt 1), 051916–051916. doi: 10.1103/PhysRevE.67.051916

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chacron, M. J., Longtin, A., St-Hilaire, M., and Maler, L. (2000). Suprathreshold stochastic firing dynamics with memory in P-type electroreceptors. Phys. Rev. Lett. 85, 1576. doi: 10.1103/PhysRevLett.85.1576

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chacron, M. J., Lindner, B., and Longtin, A. (2004). Noise shaping by interval correlations increases neuronal information transfer. Phys. Rev. Lett. 92, 080601. doi: 10.1103/PhysRevLett.92.080601

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Chacron, M. J., Longtin, A., and Maler, L. (2001). Negative interspike interval correlations increase the neuronal capacity for encoding time-dependent stimuli. J. Neurosci. 21, 5328.

Pubmed Abstract | Pubmed Full Text

Chacron, M. J., Pakdaman, K., and Longtin, A. (2003). Interspike interval correlations, memory, adaptation, and refractoriness in a leaky integrate-and-fire model with threshold fatigue. Neural Comp. 15, 253. doi: 10.1162/089976603762552915

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cox, D. R., and Lewis, P. A. W. (1966). The Statistical Analysis of Series of Events, chapter 4.6. (London: Chapman and Hall). doi: 10.1007/978-94-011-7801-3

CrossRef Full Text

Deemyad, T., Kroeger, J., and Chacron, M. J. (2012). Sub- and suprathreshold adaptation currents have opposite effects on frequency tuning. J. Physiol. 59(Pt 19), 4839–4858. doi: 10.1113/jphysiol.2012.234401

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Engel, T. A., Schimansky-Geier, L., Herz, A., Schreiber, S., and Erchova, I. (2008). Subthreshold Membrane-Potential resonances shape Spike-Train patterns in the entorhinal cortex. J. Neurophysiol. 100, 1576. doi: 10.1152/jn.01282.2007

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ermentrout, B., Pascal, M., and Gutkin, B. (2001). The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. Neural Comp. 13, 1285–1310. doi: 10.1162/08997660152002861

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ermentrout, G. B. (1996). Type I membranes, phase resetting curves, and synchrony. Neural Comp. 8, 979. doi: 10.1162/neco.1996.8.5.979

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ermentrout, G. B., and Terman, D. H. (2010). Mathematical Foundations of Neuroscience. Springer. doi: 10.1007/978-0-387-87708-2

CrossRef Full Text

Farkhooi, F., Muller, E., and Nawrot, M. P. (2011). Adaptation reduces variability of the neuronal population code. Phys. Rev. E 83(5 Pt 1), 050905–050905. doi: 10.1103/PhysRevE.83.050905

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Farkhooi, F., Strube-Bloss, M. F., and Nawrot, M. P. (2009). Serial correlation in neural spike trains: Experimental evidence, stochastic modeling, and single neuron variability. Phys. Rev. E 79, 021905–021910. doi: 10.1103/PhysRevE.79.021905

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fisch, K., Schwalger, T., Lindner, B., Herz, A., and Benda, J. (2012). Channel noise from both slow adaptation currents and fast currents is required to explain spike-response variability in a sensory neuron. J. Neurosci. 32, 17332–17344. doi: 10.1523/JNEUROSCI.6231-11.2012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Fourcaud-Trocmé, N., Hansel, D., van Vreeswijk, C., and Brunel, N. (2003). How spike generation mechanisms determine the neuronal response to fluctuating inputs. J. Neurosci. 23, 11628–11640.

Pubmed Abstract | Pubmed Full Text

Gabbiani, F., and Krapp, H. G. (2006). Spike-frequency adaptation and intrinsic properties of an identified, looming-sensitive neuron. J. Neurophysiol. 96, 2951–2962. doi: 10.1152/jn.00075.2006

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Galán, R. F., Ermentrout, G. B., and Urban, N. N. (2005). Efficient estimation of phase-resetting curves in real neurons and its significance for neural-network modeling. Phys. Rev. Lett. 94, 158101. doi: 10.1103/PhysRevLett.94.158101

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Geisler, C., and Goldberg, J. M. (1966). A stochastic model of the repetitive activity of neurons. Biophys. J. 6, 53–69. doi: 10.1016/S0006-3495(66)86639-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gerstner, W., and Naud, R. (2009). Neuroscience. how good are neuron models? Science 326, 379–380. doi: 10.1126/science.1181936

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Izhikevich, E. M. (2003). Simple model of spiking neurons. IEEE Trans. Neural. Netw. 14, 1569–1572. doi: 10.1109/TNN.2003.820440

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Izhikevich, E. M. (2005). Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. (Cambridge, MA; London: MIT Press).

Ladenbauer, J., Augustin, M., Shiau, L., and Obermayer, K. (2012). Impact of adaptation currents on synchronization of coupled exponential integrate-and-fire neurons. PLoS Comput. Biol. 8:e1002478. doi: 10.1371/journal.pcbi.1002478

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lindner, B. (2004). Interspike interval statistics of neurons driven by colored noise. Phys. Rev. E 69, 022901. doi: 10.1103/PhysRevE.69.022901

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, Y.-H., and Wang, X.-J. (2001). Spike-frequency adaptation of a generalized leaky integrate-and-fire model neuron. J. Comp. Neurosci. 10, 25. doi: 10.1023/A:1008916026143

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lowen, S. B., and Teich, M. C. (1992). Auditory-nerve action potentials form a nonrenewal point process over short as well as long time scales. J. Acoust. Soc. Am. 92, 803. doi: 10.1121/1.403950

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Middleton, J. W., Chacron, M. J., Lindner, B., and Longtin, A. (2003). Firing statistics of a neuron model driven by long-range correlated noise. Phys. Rev. E 68, 021920. doi: 10.1103/PhysRevE.68.021920

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Muller, E., Buesing, L., Schemmel, J., and Meier, K. (2007). Spike-frequency adapting neural ensembles: Beyond mean adaptation and renewal theories. Neural Comp. 19, 2958–3110. doi: 10.1162/neco.2007.19.11.2958

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Naud, R., and Gerstner, W. (2012). Coding and decoding with adapting neurons: a population approach to the peri-stimulus time histogram. PLoS Comput. Biol. 8:e1002711. doi: 10.1371/journal.pcbi.1002711

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Naud, R., Marcille, N., Clopath, C., and Gerstner, W. (2008). Firing patterns in the adaptive exponential integrate-and-fire model. Biol. Cybern. 99, 335–347. doi: 10.1007/s00422-008-0264-7

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nawrot, M. P., Boucsein, C., Rodriguez-Molina, V., Aertsen, A., Grun, S., and Rotter, S. (2007). Serial interval statistics of spontaneous activity in cortical neurons in vivo and in vitro. Neurocomputing 70, 1717. doi: 10.1016/j.neucom.2006.10.101

CrossRef Full Text

Neiman, A., and Russell, D. F. (2001). Stochastic biperiodic oscillations in the electroreceptors of paddlefish. Phys. Rev. Lett. 86, 3443. doi: 10.1103/PhysRevLett.86.3443

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Neiman, A., and Russell, D. F. (2004). Two distinct types of noisy oscillators in electroreceptors of paddlefish. J. Neurophysiol. 92, 492–509. doi: 10.1152/jn.00742.2003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Neiman, A., and Russell, D. F. (2005). Models of stochastic biperiodic oscillations and extended serial correlations in electroreceptors of paddlefish. Phys. Rev. E 71, 061915. doi: 10.1103/PhysRevE.71.061915

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nesse, W. H., Maler, L., and Longtin, A. (2010). Biophysical information representation in temporally correlated spike trains. Proc. Natl. Acad. Sci. U.S.A. 107, 21973–21978. doi: 10.1073/pnas.1008587107

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Nesse, W. H., Negro, C. A., and Bressloff, P. C. (2008). Oscillation regularity in noise-driven excitable systems with multi-time-scale adaptation. Phys. Rev. Lett. 101, 088101–088101. doi: 10.1103/PhysRevLett.101.088101

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Netoff, T., Schwemmer, M. A., and Lewis, T. J. (2012). “Experimentally estimating phase response curves of neurons: theoretical and practical issues.” in Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis, eds N. W. Schultheiss, A. A. Prinz, and R. J. Butera (Springer), 95–130.

Nikitin, A. P., Stocks, N. G., and Bulsara, A. R. (2012). Enhancing the resolution of a sensor via negative correlation: a biologically inspired approach. Phys. Rev. Lett. 109, 238103. doi: 10.1103/PhysRevLett.109.238103

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pozzorini, C., Naud, R., Mensi, S., and Gerstner, W. (2013). Temporal whitening by power-law adaptation in neocortical neurons. Nat. Neurosci. 16, 942–948. doi: 10.1038/nn.3431

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Prescott, S. A., and Sejnowski, T. J. (2008). Spike-rate coding and spike-time coding are affected oppositely by different adaptation mechanisms. J. Neurosci. 28, 13649–13661. doi: 10.1523/JNEUROSCI.1792-08.2008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ratnam, R., and Goense, J. B. M. (2004). “Variance stabilization of spike trains via non-renewal mechanisms - the impact on the speed and reliability of signal detection,” in Computational Neuroscience Meeting (CNS*2004) (Baltimore).

Ratnam, R., and Nelson, M. E. (2000). Nonrenewal statistics of electrosensory afferent spike trains: implications for the detection of weak sensory signals. J. Neurosci. 20, 6672.

Pubmed Abstract | Pubmed Full Text

Richardson, M. J. E. (2009). Dynamics of populations and networks of neurons with voltage-activated and calcium-activated currents. Phys. Rev. E 80, 021928. doi: 10.1103/PhysRevE.80.021928

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Schwalger, T., Fisch, K., Benda, J., and Lindner, B. (2010). How noisy adaptation of neurons shapes interspike interval histograms and correlations. PLoS Comput. Biol. 6:e1001026. doi: 10.1371/journal.pcbi.1001026

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Schwalger, T., and Lindner, B. (2010). Theory for serial correlations of interevent intervals. Eur. Phys. J. Spec. Top. 187, 211–221. doi: 10.1140/epjst/e2010-01286-y

CrossRef Full Text

Softky, W. R., and Koch, C. (1993). The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J. Neurosci. 13, 334.

Pubmed Abstract | Pubmed Full Text

Treves, A. (1993). Mean-field analysis of neuronal spike dynamics. Netw. Comput. Neural Syst. 4, 259–284. doi: 10.1088/0954-898X/4/3/002

CrossRef Full Text

Urdapilleta, E. (2011). Onset of negative interspike interval correlations in adapting neurons. Phys. Rev. E 84, 041904. doi: 10.1103/PhysRevE.84.041904

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van Vreeswijk, C. (2010). “Stochastic models of spike trains,” in Analysis of Parallel Spike Trains, eds S. Grün and S. Rotter (Springer).

Vogel, A., Hennig, R. M., and Ronacher, B. (2005). Increase of neuronal response variability at higher processing levels as revealed by simultaneous recordings. J. Neurophysiol. 93, 3548. doi: 10.1152/jn.01288.2004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wang, X. J. (1998). Calcium coding and adaptive temporal computation in cortical pyramidal neurons. J. Neurophysiol. 79, 1549–1566.

Pubmed Abstract | Pubmed Full Text