Abstract

The property of a neuron to phase-lock to an oscillatory stimulus before adapting its spike rate to the stimulus frequency plays an important role for the auditory system. We investigate under which conditions neurons exhibit this phase locking below rate threshold. To this end, we simulate neurons employing the widely used leaky integrate-and-fire (LIF) model. Tuning parameters, we can arrange either an irregular spontaneous or a tonic spiking mode. When the neuron is stimulated in both modes, a significant rise of vector strength prior to a noticeable change of the spike rate can be observed. Combining analytic reasoning with numerical simulations, we trace this observation back to a modulation of interspike intervals, which itself requires spikes to be only loosely coupled. We test the limits of this conception by simulating an LIF model with threshold fatigue, which generates pronounced anticorrelations between subsequent interspike intervals. In addition we evaluate the LIF response for harmonic stimuli of various frequencies and discuss the extension to more complex stimuli. It seems that phase locking below rate threshold occurs generically for all zero mean stimuli. Finally, we discuss our findings in the context of stimulus detection.

1.  Introduction

Phase locking—the firing of spikes at a preferred phase of an oscillatory stimulus—is a universal mechanism in sensory neurons and can be observed for nearly all sensory modalities. Early reports exist for cat retinal ganglion cells (Ogawa, Bishop, & Levick, 1966), single auditory nerve fibers of the squirrel monkey (Rose, Brugge, Anderson, & Hind, 1967), and mechanoreceptive afferents of the monkey hand (Talbot, Darian-Smith, Kornhuber, & Mountcastle, 1968). In the olfactory system, a phase locking of mitral cells to intrinsic subthreshold oscillations is discussed as an essential mechanism for robust odor representation (Desmaisons, Vincent, & Lledo, 1999; Wilson & Mainen, 2006). Finally, different types of noisy phase-locked regimes were observed in the electroreceptor of the paddlefish stimulated by oscillating electric dipole fields (Neiman et al., 1999).

Phase locking and synchronization are fundamental concepts playing an important role in neuronal information processing systems. Although both effects are related (phase locking is sometimes termed synchronization in the statistical sense), it is important to clearly distinguish between them: whereas phase locking means that in a drive response relation, the response occurs at a preferential phase of the driving signal, synchronization (Pikovsky, Rosenblum, & Kurths, 2001) is characterized by the forced matching of spike rate and stimulus frequency or the boundedness of a suitably defined phase difference. Eliciting a spike at each peak of a sinusoidal stimulus describes a situation where phase locking and synchronization co-occur. Spikes riding on the signal peaks but skipping cycles preserve phase locking but introduce phase slips, thus corrupting synchronization. On the other hand, randomizing the position of each spike in its respective periodic interval would erase phase locking but preserve synchronization.

In case of a stimulus with vanishing strength, neither phase locking nor synchronization will be observed. Gradually increasing the stimulus strength often leads to the observation that significant phase locking sets in before a noticeable change of the spontaneous spike rate, that is, before reaching rate threshold. This phase locking below rate threshold is particularly prominent in the auditory nerve (Johnson, 1980; Köppl, 1997). It was also found in a conceptual model aiming to explain the exquisite sensitivity of the auditory system to acoustic stimuli that carry less energy than the background noise (Camalet, Duke, Jülicher, & Prost, 2000, Fig. 6). This raises the question of what properties of a neuron are necessary or sufficient to effect phase locking below rate threshold. Surprisingly, we did not find an answer to this very basic question in the existing literature. Here we try to close the gap through a modeling approach and dwell on the mechanism underlying this effect.

2.  Methods

We combine two methods to elucidate the conditions under which phase locking below rate threshold occurs: (1) we derive a mathematical relation between the spike rate and the vector strength, the latter being a common measure for phase locking; and (2), we perform diverse numerical simulations of the driven leaky-integrate-and-fire neuron model.

2.1.  Relation Between Spike Rate and Vector Strength.

We consider a harmonic stimulus s(t) = a cos(ϕ0 + Ωt) with a linearly increasing phase ϕ(t) = ϕ0 + Ωt, frequency Ω and strength a. The restriction to zero mean harmonic stimuli is motivated by the filter bank property of the cochlea (Glasberg & Moore, 1990; Hohmann, 2002) and the fact that phase locking below rate threshold is most prominent in the auditory system. The case of more complex stimuli is addressed in section 3.2.

The zero mean harmonic stimulus drives a model neuron that, according to some specific rules, elicits spikes at times ti. Differences between these spike times constitute the sequence of interspike intervals (ISIs) τi = titi−1. Moreover, we collect the wrapped stimulus phases at the spike times ϕi = ϕ(ti)(mod 2π). Binning these spike phases, we construct the so-called period histogram and, after normalization and division by the bin widths, estimate a related probability density p(ϕ). This phase distribution is a central quantity for a detection of phase locking seen in a concentration at a certain value. The degree of phase locking can be quantified by the so-called vector strength (or synchronization index) (Goldberg & Brown, 1969):
formula
2.1
This quantity varies between zero for the equidistribution (if p(ϕ) = (2π)−1), and one for a delta distribution (if p(ϕ) = δ(ϕ − ϕ0)). Using a Fourier decomposition of the 2π-periodic function p(ϕ),
formula
2.2
with Fourier coefficients,
formula
2.3
we see immediately that the vector strength is nothing but the modulus of the first Fourier coefficient: v = |c1| = |c−1|. Notice that c0 = 1 due to normalization of p(ϕ) and ck = c*k since p(ϕ) is real. Moreover, we note that via the phase distribution p(ϕ) = p(ϕ; a, Ω), all Fourier coefficients are functions of the stimulus strength and frequency (i.e., ck = ck(a, Ω)) and, in particular, v = v(a, Ω).
The second quantity we focus on is the spike rate r of the driven neuron. Here, the term spike rate relates to spike counts, not to an ensemble-based instantaneous firing rate (or peristimulus time histogram PSTH).1 Its operative definition requires counting the number N of spikes a neuron fires within the time window T and then computing the ratio
formula
2.4
Taking the inverse of this relation and noticing that the time window T can be chosen to encompass (N − 1) consecutive ISIs τ1, τ2, …, τN−1, we find the relation
formula
2.5
Here we have assumed that downstream neurons compute the inverse spike rate by measuring the time passing between an initial and an Nth spike. In case the downstream neuron estimates by counting the number of spikes observed in a preselected fixed time window T, equation 2.5, would require a correction factor where 〈τ〉 denotes the average ISI.
For N → ∞ the expression in equation 2.5 can be interpreted as a mean ISI 〈τ〉 where the long-time statistics is described by a probability density p(τ; a, Ω) that covaries with the stimulus parameters a and Ω. Hence, we find
formula
2.6
In principle we should also account for the initial phase ϕ0 of the stimulus and perform an additional average. However, in the neurons we have in mind, the memory for this initial phase decays and is lost beyond a transient regime (ergodicity assumption).
We link the quantity 2.6 to the vector strength via the identity
formula
2.7
which leads to
formula
2.8
formula
2.9
formula
2.10
and by insertion of the Fourier series, equation 2.2,
formula
2.11
formula
2.12
formula
2.13
where we have employed the definitions
formula
2.14
formula
2.15
and the notation ℜ{·} for the real part of a complex number. Writing and (and suppressing further notation of stimulus parameters a and Ω), we can rewrite equation 2.13 as
formula
2.16
For vanishing stimulus strength (a = 0), the phase distribution p(ϕ) is flat, corresponding to ck⩾1 = 0, and ϑ0 = 〈τ〉. Increasing the stimulus strength causes a structuring of the phase distribution that is reflected in the Fourier coefficients ck>0 gaining increasing weight. The spectrum of Fourier coefficients exhibits a decaying profile: the larger k, the smaller the modulus of ck. The background of this observation is the reciprocal relation between the effective widths in original ϕ-space and its related Fourier space, well known from quantum mechanical wave packets and the uncertainty relation. Notice that the singular zero-width delta distribution p(ϕ) = δ(ϕ − ϕ0) corresponds to the infinite-width spectrum ck = 1 for all k. In line with this observation, a weak stimulus (ws) approximation is expressed as (ck⩾2 ≈ 0) and means that the initial deviation of p(ϕ) from the flat distribution is nothing but a harmonic modulation with the period 2π, that is,
formula
2.17
and with the peak located at the phase γ1. Because of the nonnegativity of the phase distribution, formula 2.17 must necessarily be restricted to values v ⩽ 0.5. Moreover, for weak stimuli the relation between spike rate r and vector strength v, equation 2.16, can be approximated as
formula
2.18
Now we can see how an increase in vector strength, say up to a value of 0.5, can occur without being accompanied by a change in spike rate. This phase locking below rate threshold stems from the fact that ϑ0 remains effectively unchanged and |ϑ1| ≪ ϑ0. Equation 2.16 shows that extending the second property to |ϑk| ≪ ϑ0 for all k ⩾ 1 even yields a sufficient criterion for phase locking below rate threshold without restriction to weak stimuli.

2.2.  Noisy Leaky-Integrate-and-Fire Neuron.

One of the simplest and most widely used models for neuronal spike generation is the leaky-integrate-and-fire (LIF) model (Stein, 1965). While the integrate-and-fire model is usually ascribed to Lapicque and his often cited 1907 paper (Lapicque, 1907; Brunel & van Rossum, 2007b), it has recently been pointed out (Brunel & van Rossum, 2007a) that it was in fact introduced by theorists much later.

The LIF model is based on splitting a synaptic current Is into a resistive leak current and a capacitive current , where u(t), R, and C are membrane voltage, resistor, and capacitance, respectively. Rearranging the equation Is = IR + IC, we can write
formula
2.19
with the membrane time constant τm = RC. One component of the synaptic current results from a bombardment of the LIF neuron with spikes arriving from a multitude of presynaptic neurons. Under certain assumptions (Lansky & Sato, 1999) its effect can be modeled by a mean current I0 plus gaussian white noise ζ(t). Differing from this frequently considered case, gaussian red (or correlated) noise accounts for synchronous presynaptic firing; its effects on the ISI statistics and serial correlations of a LIF model were recently investigated by Schwalger and Schimansky-Geier (2008). As a second component of the synaptic current, we add a harmonic stimulus with strength , frequency Ω, and initial phase ϕ0 = π (cf. the ergodicity assumption above). Insertion in equation 2.19yields
formula
2.20
Next we add a spike triggering voltage threshold uθ and a reset voltage ur. We complete our LIF model neuron by endowing it with an absolute refractory time τa = 0.5 during which integration ceases. Choosing ur as the voltage reference and going to a rescaled membrane voltage x(t) = [u(t) − ur]/uθ, we finally arrive at
formula
2.21
Here the first term describes the deterministic relaxation: in the absence of second (σ = 0) and third (a = 0) terms, the voltage variable x(t) relaxes to the value x within a time of the order of a few τm. The second term in equation 2.21 reflects rescaled voltage fluctuations in the form of standard zero mean additive gaussian white noise,
formula
2.22
of intensity σ2. The first two terms in equation 2.21 make up what is called an Ornstein-Uhlenbeck process (Gardiner, 2004). The third term in equation 2.21 expresses the rescaled harmonic stimulus.

By choice of x, we can arrange two different operational modes of our model neuron. Placing the resting voltage x slightly below the firing threshold xθ (case A) leaves a gap xθx that has to be traversed with the aid of random fluctuations. This situation corresponds to noise-mediated spontaneous spiking activity. On the other hand, fixing x above xθ (case B) leads to tonic spiking, with the width of the ISI distribution being controlled by the noise intensity σ.

We point out that by construction as a renewal process, serial ISI correlations in the standard LIF are absent. A physiologically plausible way to implement these experimentally observed negative ISI correlations is to reformulate the LIF with threshold fatigue (Chacron, Pakdaman, & Longtin, 2003). This means the transition to a dynamic threshold θ(t), which is relaxing with time constant τθ toward xθ (case C):
formula
2.23
As before, the dynamics of x(t) is governed by equation 2.21. The very moment that x(t) hits θ(t) from below, x(t) is reset to xr whereas θ(t) is set to θ0 + W[θ(t), α]. While in Chacron et al. (2003), W was required to be a monotonically increasing function of both arguments obeying W[θ(t), 0] = 0, here we will consider only the special case W[θ(t), α = 1] = θ(t) and θ0 = 1. This means that each time a spike is fired, the threshold is incremented by unity. We note that the standard LIF with constant threshold xθ = 1 corresponds to θo = xθ and α = 0.

All three cases A, [B], and {C}, were simulated numerically employing the Heun scheme (Mannella, 2000). The results, depicted in Figure 1, were obtained with parameters given in the figure caption, chosen to yield a similar spontaneous spike rate for all three cases.

Figure 1:

Three operating modes of the undriven noisy leaky-integrate-and-fire model. (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking induced via a dynamic threshold. Parameters relating to equations 2.21 and 2.23 are chosen to yield comparable spike rates r ≈ 0.12 in all three cases A[B]{C}: reset value xr = 0, static threshold xθ = 1, asymptote x = 0.9[1.1]{2}, absolute refractory period τa = 0.5, membrane relaxation time τm = 0.5[3.3]{0.5}, noise strength σ = 0.1[0.025]{0.5}, threshold relaxation time {τθ = 15}.

Figure 1:

Three operating modes of the undriven noisy leaky-integrate-and-fire model. (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking induced via a dynamic threshold. Parameters relating to equations 2.21 and 2.23 are chosen to yield comparable spike rates r ≈ 0.12 in all three cases A[B]{C}: reset value xr = 0, static threshold xθ = 1, asymptote x = 0.9[1.1]{2}, absolute refractory period τa = 0.5, membrane relaxation time τm = 0.5[3.3]{0.5}, noise strength σ = 0.1[0.025]{0.5}, threshold relaxation time {τθ = 15}.

3.  Results

3.1.  Numerical Simulation Results of the Driven LIF Model Neuron.

In the following we present simulation results of the driven noisy LIF for the stimulus frequency Ω = 2π and amplitudes ranging from a = 0.001, …, 10 on a decadal logarithmic scale. The other parameters for the cases A[B]{C} are as given in the caption of Figure 1. The mean spike rate r0 in the absence of the stimulus is roughly 0.12 for all cases; hence, the mean ISI is roughly 8.3 times larger than the period of the stimulus Ts = 2π/Ω = 1. This means that the transition from the undriven neuronal dynamics via phase locking to the (1:1) synchronization (one spike per stimulus period) regime is accompanied by an uprise of the spike rate—to synchronize, the neuron has to speed up. The stimulus-induced increase of vector strength v and spike rate r, computed for all three cases from numerical simulations of 20,000 spikes per amplitude, are displayed in the three panels of Figure 2, respectively.

Figure 2:

Phase locking below rate threshold in the noisy LIF can be seen in the range of amplitudes where the vector strength (black solid curve) is already elevated, while the spike rate (gray solid curve) is still near the level of spontaneous activity (r0 ≈ 0.12). The black dashed line marks the rate threshold at defined as in equation 3.1. (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking induced via a dynamic threshold. The serial ISI correlation coefficient ρ1 ≈ −0.5 shown in the inset of the bottom panel reflects a pronounced anticorrelation between subsequent ISIs, which, however, does not prevent phase locking but hampers an increase of the spike rate. Parameters are as given in the caption of Figure 1 and Ω = 2π.

Figure 2:

Phase locking below rate threshold in the noisy LIF can be seen in the range of amplitudes where the vector strength (black solid curve) is already elevated, while the spike rate (gray solid curve) is still near the level of spontaneous activity (r0 ≈ 0.12). The black dashed line marks the rate threshold at defined as in equation 3.1. (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking induced via a dynamic threshold. The serial ISI correlation coefficient ρ1 ≈ −0.5 shown in the inset of the bottom panel reflects a pronounced anticorrelation between subsequent ISIs, which, however, does not prevent phase locking but hampers an increase of the spike rate. Parameters are as given in the caption of Figure 1 and Ω = 2π.

In the simulated noisy LIF, phase locking below rate threshold can be seen in the range of amplitudes where the vector strength (black solid curve) is already elevated, while the spike rate (gray solid curve) is still near the level of spontaneous activity (r0 ≈ 0.12). Quantifying the rate threshold as
formula
3.1
and the degree of phase locking at rate threshold by v(at), we find for case A, at ≈ 0.1, v(at) ≈ 0.4; for case B, at ≈ 0.25, v(at) ≈ 0.85; and for case C: at ≈ 2, v(at) ≈ 0.81. In the viewgraphs of Figure 2, these items are indicated by the black dashed line.

To gain further insight into the phenomenon and the structure of the curves, we also show plots of the ISI and phase histograms in Figure 3. Notice that the gray lines trace the mean ISI (or inverse spike rate) and the vector strength, respectively.

Figure 3:

(Left column): Gray-coded ISI histograms (along abscissa) for a wide range of stimulus amplitudes (along ordinate). Gray solid lines trace the mean ISI or the inverse spike rate. (Right column): Gray-coded phase histograms (along abscissa) for a wide range of stimulus amplitudes (along ordinate). Gray solid lines trace the vector strength (cf. Figure 2). (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking. Parameters are as given in the caption of Figure 1 and Ω = 2π.

Figure 3:

(Left column): Gray-coded ISI histograms (along abscissa) for a wide range of stimulus amplitudes (along ordinate). Gray solid lines trace the mean ISI or the inverse spike rate. (Right column): Gray-coded phase histograms (along abscissa) for a wide range of stimulus amplitudes (along ordinate). Gray solid lines trace the vector strength (cf. Figure 2). (Top) Case A: Noise-activated spiking. (Middle) Case B: Tonic spiking. (Bottom) Case C: Correlated spiking. Parameters are as given in the caption of Figure 1 and Ω = 2π.

3.2.  More Complex Stimuli.

Studying phase locking to more complex stimuli first raises the question of how to generalize the phase concept for nonharmonic or even nonperiodic stimuli. As discussed in the context of phase synchronization (Freund, Schimansky-Geier, & Hänggi, 2003), several phase definitions (natural, linear interpolating, Hilbert phase) may qualify. After an appropriate phase definition ϕ(t) has been chosen, a general stimulus may be written as s(t) = 〈s〉 + a(t)cos ϕ(t). The harmonic case considered above corresponds to 〈s〉 = 0, a(t) = a, ϕ(t) = Ωt.

The case demonstrates two things: first, that the appropriate phase, here ϕ(t) = 2Ωt, can deviate from the naive expression Ωt suggested by the definition and, second, that nonlinearities may generate a nonvanishing mean, here 〈s〉 = a/2, which covaries with the amplitude of the harmonic oscillatory component. A nonvanishing mean breaks the balance between depolarizing and hyperpolarizing stimulus segments; if ∂|〈s〉|/∂a>0, the stronger the stimulus. This balance was underlying the pure shifting of spikes observed for the unbiased harmonic stimuli considered above. From this reasoning, it is straightforward to anticipate that the harmonic but biased stimulus changes the spike rate before noticeable phase locking sets in. Through numerical simulations shown in the top panel of Figure 4 we found that for this stimulus, phase locking below rate threshold cannot be observed.

Figure 4:

Simulating the LIF response to three stimuli sketched in the panel insets (for a = 1), we investigated whether phase locking below rate threshold occurs also for more complex stimuli: a nonzero mean harmonic stimulus (top panel) and two strongly asymmetric but zero mean stimuli (cf. equation 3.2): middle panel and bottom panel). From our simulation results, we conclude that phase locking below rate threshold can also occur for anharmonic stimuli as long as they have a vanishing mean. In contrast, the strength-dependent mean a/2 (realized in the top panel) significantly lowers the rate threshold, thus preventing phase locking below rate threshold. The multipeaked structure of the vector strength seen in the top panel is caused by recurring peak splitting. LIF parameters are as given in the caption of Figure 1 for case A.

Figure 4:

Simulating the LIF response to three stimuli sketched in the panel insets (for a = 1), we investigated whether phase locking below rate threshold occurs also for more complex stimuli: a nonzero mean harmonic stimulus (top panel) and two strongly asymmetric but zero mean stimuli (cf. equation 3.2): middle panel and bottom panel). From our simulation results, we conclude that phase locking below rate threshold can also occur for anharmonic stimuli as long as they have a vanishing mean. In contrast, the strength-dependent mean a/2 (realized in the top panel) significantly lowers the rate threshold, thus preventing phase locking below rate threshold. The multipeaked structure of the vector strength seen in the top panel is caused by recurring peak splitting. LIF parameters are as given in the caption of Figure 1 for case A.

Perhaps more interesting than a biased harmonic (or symmetric about nonzero mean) stimulus is the effect of an unbiased anharmonic (or asymmetric about zero mean) stimulus. The vanishing mean expresses a balance of the first moment but still leaves room for an imbalance of higher odd moments. As a test case, we simulated the LIF neuron driven by a strongly asymmetric stimulus,
formula
3.2
which, by construction, has a vanishing mean and oscillates between ∓0.15a and ±a. The strong anharmonicity of in the time domain corresponds to a significant weight placed on the higher harmonics in the spectral domain. Defining the instantaneous phase by linear interpolation between the regular maxima amounts to ϕ(t) = 2πt. Simulation results with and , shown in the middle and bottom panels of Figure 4, respectively reveal that in spite of their pronounced asymmetry, these zero mean stimuli nonetheless induce phase locking below rate threshold.

3.3.  Evaluation of the Relation Between Spike Rate and Vector Strength.

An evaluation of equation 2.16 or its weak stimulus approximation, equation 2.17, requires a computation of the Fourier coefficients ck of the phase distribution 2.3 and of ϑk, the Fourier coefficients of the conditional mean ISI, which occur in the expansion,
formula
3.3
and are defined in equation 2.14. While the ϑk are related to p(τ ∣ ϕ; a, Ω) via equation 2.15, the Fourier coefficients ck can be connected to the same conditional ISI probability density via Fourier coefficients of higher conditional moments ϑnk. This connection is not explicated here for the reason that a closed analytic expression for the basic ingredient p(τ ∣ ϕ; a, Ω) is, to the best of our knowledge, not available for the driven noisy LIF. Results obtained by Burkitt and Clark (2001) and Kempter, Gerstner, van Hemmen, and Wagner (1998) in closely related investigations approach desired expressions, but finally these also rely on numerical evaluation. A detailed elaboration may be feasible for the perfect integrator with random threshold (Gabbiani & Koch, 1996); this will be left to future research.
Instead, we obtained p(τ ∣ ϕ; a, Ω) through numerical simulations, the results of which are displayed in Figure 5. A pronounced structure in the (τ, ϕ)-plane (with period Ts along the τ-axis) emerges for intermediate stimulus strength, that is, at rate threshold (middle panel). This structure is reflected in the Fourier coefficients ck and, in particular, in a significant vector strength v = |c1|. To elaborate the connection, one has to express the stationary phase distribution p(ϕ; a, Ω) as the eigenfunction to the eigenvalue 1 of the transition kernel (Plesser & Geisel, 1999):
formula
3.4

The coefficients ϑk are also illustrated in Figure 5. The black and gray lines show projections of the conditional mean ISI 〈τ ∣ ϕ; a, Ω〉 for the driven and undriven LIF model, respectively. For a = 0 (bottom panel), both coincide, and, as expected, no dependence on ϕ is visible. The position of the gray line is thus given by ϑ−10(a = 0, Ω). As a approaches the rate threshold, the black line deviates from the gray line within the bound formulated by equation 3.1. From the middle panel of Figure 5, we see that the black line still exhibits no dependence on ϕ. This is the graphical equivalent of the fact that |ϑk⩾1| ≪ ϑ0 and that the black line is solely determined by ϑ−10(a, Ω). This separation of scales remains valid for intermediate values of stimulus strength where ϑ0 starts to decline; even beyond rate threshold, we find that up to a good approximation, r(a, Ω) ≈ ϑ−10(a, Ω). Only for very strong stimuli (top panel of Figure 5) do we see a pronounced modulation of the black line corresponding to a nonnegligible Fourier coefficient ϑ1(a, Ω). However, this is far beyond rate threshold in the regime of frequency and phase synchronization and not plain phase locking.

Figure 5:

The conditional ISI probability density p(τ ∣ ϕ; a, Ω) computed from numerical simulations of the driven noisy LIF model neuron (case A: all parameters except a as specified in the caption of Figure 1 and Ω = 2π). Panels bottom to top for a = 0 (undriven neuron), a = 0.1 (at rate threshold) and a = 1 (strong stimulus), respectively. Condition ϕ denotes the stimulus phase at the moment of the previous spike. The gray and black solid lines in the projection plane trace the conditional mean ISI 〈τ ∣ ϕ; a, Ω〉 for the undriven and driven neuron, respectively.

Figure 5:

The conditional ISI probability density p(τ ∣ ϕ; a, Ω) computed from numerical simulations of the driven noisy LIF model neuron (case A: all parameters except a as specified in the caption of Figure 1 and Ω = 2π). Panels bottom to top for a = 0 (undriven neuron), a = 0.1 (at rate threshold) and a = 1 (strong stimulus), respectively. Condition ϕ denotes the stimulus phase at the moment of the previous spike. The gray and black solid lines in the projection plane trace the conditional mean ISI 〈τ ∣ ϕ; a, Ω〉 for the undriven and driven neuron, respectively.

4.  Discussion

4.1.  Phase Locking and ISI Correlations.

The essence of phase locking below rate threshold can be traced back to the fact that while the probability density p(τ ∣ ϕ; aΩ) is a highly structured function of ϕ the average 〈τ ∣ ϕ; a, Ω〉 (see equation 2.15) is nearly independent from ϕ. This lack of sensitivity also explains why a substitution of the structured p(ϕ; a, Ω) in equation 2.9 by the equidistribution [2π]−1 does not make a big difference, thus justifying the approximation r(a, Ω) ≈ ϑ−10(a, Ω) valid for weak to intermediate stimuli (to see this connection, specify equation 2.14 for k = 0).

The structure in p(τ ∣ ϕ; aΩ) and the independence of 〈τ ∣ ϕ; a, Ω〉 from ϕ are caused by a precise balance of spikes advanced or delayed by the stimulus to match a favorable stimulus phase. Phase locking below rate threshold thus means eliciting the same average number of spikes that would occur in the absence of the stimulus and establishing this balance by shifting (i.e., advancing or delaying) these spikes. In synchronization research, this mechanism is also known as modulation and contrasted with synchronization (Pikovsky et al., 2001). In particular, this means that subsequent spikes can be shifted independently, which raises the question of what kind of correlations prevent phase locking below rate threshold.

To see how interspike correlations can affect the phase-locking properties of a neuron, we imagine all spikes to be stiffly coupled, like the teeth of a comb, with a strong tendency to maintain fixed interspike intervals. Shifting the first spike to match a preferred phase generically would lead to a mismatch of some of the following spikes. Phase locking—the excess match of spikes with a preferred stimulus phase—could then be arranged only by excess insertion (or excess deletion) of spikes, that is, a change of the spike rate. Notice that pure spike shifting to the preferred phase (i.e., phase locking below rate threshold) can always be viewed as a balanced deletion (at the mismatch phase) and insertion (at the preferred phase).

Correlations between subsequent spikes are usually discussed in the context of so-called serial ISI correlations that have been observed in spike trains of sensory neurons (e.g., Bahar et al., 2001). A prominent anticorrelation between subsequent ISIs—a short ISI is followed by a long one and vice versa—was reported for the P-type electroreceptor of the weakly electric fish (Ratnam & Nelson, 2000; Chacron, Longtin, St.-Hilaire, & Maler, 2000). Suitable models reproducing these effects are the LIF with threshold fatigue (Chacron et al., 2003) or a perfect integrator with a noisy threshold and a reset by a constant decrement (Chacron, Lindner, & Longtin, 2004; Lindner, Chacron, & Longtin, 2005). Simulation results of a LIF with threshold fatigue are depicted in the bottom panel of Figure 2. The serial correlation coefficient ρ1 ≈ −0.5 in the inset clearly expresses a pronounced anticorrelation. Surprisingly, this strong coupling between subsequent ISIs does not prevent phase locking below rate threshold but, on the contrary, impedes an increase of the spike rate.

4.2.  Dependence on Stimulus Frequency.

Phase locking and frequency synchronization in our periodically driven model neuron are based on an interplay of diverse timescales: the intrinsic absolute refractory time τa; the intrinsic relaxation time τm, which together and in conjunction with the noise intensity σ2 and the gap xθx control the average ISI; in the LIF with threshold fatigue the intrinsic relaxation time τs of the threshold; and, finally, the external stimulus period TΩ = 2π/Ω. Keeping all intrinsic timescales fixed, we varied the stimulus frequency Ω and computed the rate threshold at (see equation 3.1) and the vector strength at rate threshold v(at). Simulation results for both cases A and B are depicted in Figure 6. To interpret the spectral profiles, we note the following facts: there is a general trend of a rising rate threshold at with increasing stimulus frequency Ω, meaning that fast-oscillating stimuli do not change spike counts so easily. While the tonic spiking mode (case B) displays resonant structures—elevated rate thresholds at ratios Ω:2πr0 = (1:1), (1:2), (2:1)—the noise-activated spiking mode (case A) exhibits barely any structure (besides the aforementioned mentioned trend). The vector strength at rate threshold is, on the one hand, correlated with the rate threshold curve: the higher at, the higher v(at). On the other hand, for strong stimuli, there is a tendency of declining vector strength with increasing stimulus frequency. The reason for this decline is found in the noise-induced temporal jitter Δt of spikes, which translates into a dispersion of related phases via Δϕ = ΩΔt. Following Hill, Stange, and Mo (1989), one can derive a parabola v(Ω) = 1 − (ΩΔt)2/2 that approximately describes this high-frequency decline of the vector strength (up to ). For the low-frequency branch, there is a reduction of the vector strength due to peak splitting of the phase histogram (Köppl, 1997) caused by more than just one spike being fired at the slowly changing preferred stimulus phase.

Figure 6:

The dependence of vector strength at rate threshold on the frequency of the stimulus for (case A) noise-activated spiking: circles and full line, and for (case B) tonic spiking: diamonds and dashed line.

Figure 6:

The dependence of vector strength at rate threshold on the frequency of the stimulus for (case A) noise-activated spiking: circles and full line, and for (case B) tonic spiking: diamonds and dashed line.

4.3.  Phase Locking as a Detection Mechanism.

While a rate threshold at and the related vector strength v(at) are clearly defined quantities, the question of what to call phase locked and what not is not so easy to answer. A value of, say, 0.4 is not specific for a narrow single-peaked phase distribution and ambiguously may result from a variety of differently shaped distributions. Nevertheless, this value still might be sufficient for a task like sound localization. All that matters is whether the computed empirical vector strength,
formula
4.1
is safely above values that can be expected from the statistics in the absence of a stimulus (null hypothesis). Under the assumption of N independent and identically equidistributed phases ϕ1, …, ϕN (i.e., for undriven neurons modeled by a renewal process), this statistics is described by the Rayleigh distribution . The probability of observing a value of is thus given by exp(−N 0.16) (Naundorff & Freund, 2002), which crosses the 5% (1%) significance level after N = 20(30) spikes, respectively. Of course, still larger values, say, , allow rejection of the null hypothesis of stimulus absence after 5(7) spikes, respectively.2

5.  Summary and Outlook

The starting point of our investigation was the question of what specific ingredients of a neuron (model) are necessary to observe phase locking to an oscillatory stimulus before significant changes of the spike rate. Through numerical simulations, we found that even in a very basic neuron model, the noisy LIF, phase locking below rate threshold actually does occur without the need to tune parameters. Moreover, we found that it is not tied to a certain spiking mode but happens generically for both noise-induced spiking and tonic spiking. From analytic considerations, we were able to trace phase locking below rate threshold back to the neuron's inherent freedom to shift subsequent spikes independently. This insight led us to conjecture that pronounced interspike correlations might be a mechanism to preclude phase locking below rate threshold. Surprisingly, we found that significant anticorrelations between subsequent interspike intervals (ρ1 ≈ −0.5), as they occur in the LIF with threshold fatigue, are not enough to prevent phase locking but rather stiffen the spontaneous spiking rate. Also, a variation of the stimulus frequency never led to a drastic collapse of the vector strength at rate threshold, as one might expect for driving the tonic spiking neuron at half the spontaneous rate. Paradoxically it seems that specific ingredients of a neuron model are required not to guarantee but, on the contrary, to prevent phase locking to zero mean stimuli below rate threshold. Finally, we discussed the ambiguous relation between phase distributions and moderate values of the vector strength that still might qualify for a detection scheme.

In order to fully exploit the analytic relation between vector strength and spike rate explicated in this letter, it seems desirable to find a neuron model that allows deriving a closed-form expression for the core piece p(τ|ϕ; a, Ω). A good candidate seems to be the perfect integrator with stochastic threshold and reset; however, detailed elaboration is left to future research. Eventually further studies are needed to construct a reasonably realistic neuron model where phase locking to zero mean stimuli occurs only in combination with a significant variation of the spike rate.

Acknowledgments

J.A.F. acknowledges support by the DAAD (ARC D/06/12701). N.G.S. acknowledges support by grant EP/DO51894. We enjoyed fruitful discussions with S. Beulig, M. Buschermöhle, U. Feudel, R. P. Morse, and L. Schimansky-Geier.

Note

1

An account of the differences is given in Rieke, Warland, de Ruyter van Steveninck, and Bialek (1997).

2

We note that the Rayleigh distribution results from the central limit theorem; hence, for too few spikes, the statement has to be taken cum grano salis.

References

Bahar
,
S.
,
Kantelhardt
,
J. W.
,
Neiman
,
A.
,
Rego
,
H. H. A.
,
Russell
,
D. F.
,
Wilkens
,
L.
, et al
(
2001
).
Long-range temporal anti-correlations in paddlefish electroreceptors
.
Europhys. Lett.
,
56
,
454
460
.
Brunel
,
N.
, &
van Rossum
,
M. C. W.
(
2007a
).
Lapicque's 1907 paper: From frogs to integrate-and-fire
.
Biol. Cybern.
,
97
,
337
339
.
Brunel
,
N.
, &
van Rossum
,
M. C. W.
(
2007b
).
Quantitative investigations of electrical nerve excitation treated as polarization
.
Biol. Cybern.
,
97
,
341
349
.
Burkitt
,
A. N.
, &
Clark
,
G. M.
(
2001
).
Synchronization of the neural response to noisy periodic synaptic input
.
Neural Comput.
,
13
,
2639
2672
.
Camalet
,
S.
,
Duke
,
T.
,
Jülicher
,
F.
, &
Prost
,
J.
(
2000
).
Auditory sensitivity provided by self-tuned critical oscillations of hair cells
.
PNAS
,
97
,
3183
3188
.
Chacron
,
M. J.
,
Lindner
,
B.
, &
Longtin
,
A.
(
2004
).
Noise shaping by interval correlations increases information transfer
.
Phys. Rev. Lett.
,
92
,
080601
1-4
.
Chacron
,
M. J.
,
Longtin
,
A.
,
St. Hilaire
,
M.
, &
Maler
,
L.
(
2000
).
Suprathreshold stochastic firing dynamics with memory in P-type electroreceptors
.
Phys. Rev. Lett.
,
85
,
1576
1579
.
Chacron
,
M. J.
,
Pakdaman
,
K.
, &
Longtin
,
A.
(
2003
).
Interspike interval correlations, memory, adaptation, and refractoriness in a leaky integrate-and-fire model with threshold fatigue
.
Neural Comput.
,
15
,
253
278
.
Desmaisons
,
D.
,
Vincent
,
J.-D.
, &
Lledo
,
P.-M.
(
1999
).
Control of action potential timing by intrinsic subthreshold oscillations in olfactory bulb output neurons
.
J. Neurosci.
,
19
,
10727
10737
.
Freund
,
J. A.
,
Schimansky-Geier
,
L.
, &
Hänggi
,
P.
(
2003
).
Frequency and phase synchronization in stochastic systems
.
Chaos
,
13
,
225
238
.
Gabbiani
,
F.
, &
Koch
,
C.
(
1996
).
Coding of time-varying signals in spike trains of integrate-and-fire neurons with random threshold
.
Neural Comput.
,
8
,
44
66
.
Gardiner
,
C. W.
(
2004
).
Handbook of Stochastic Methods
(3rd ed.).
Berlin
:
Springer
.
Glasberg
,
B. R.
, &
Moore
,
B. C. J.
(
1990
).
Derivation of auditory filter shapes from notched-noise data
.
Hear. Res.
,
47
,
103
138
.
Goldberg
,
J. M.
, &
Brown
,
P. B.
(
1969
).
Response of binaural neurons of dog superior olivary complex to dichotic tonal stimuli: Some physiological mechanisms of sound localization
.
J. Neurophysiol
,
32
,
613
636
.
Hill
,
K. G.
,
Stange
,
G.
, &
Mo
,
J.
(
1989
).
Temporal synchronization in the primary auditory response in the pigeon
.
Hear. Res.
,
39
,
63
74
.
Hohmann
,
V.
(
2002
).
Frequency analysis and synthesis using a gammatone filterbank
.
Acta Acust. United Ac.
,
88
,
433
442
.
Johnson
,
D. H.
(
1980
).
The relationship between spike rate and synchrony in responses of auditory-nerve fibers to single tones
.
J. Acoust. Soc. Am.
,
68
,
1115
1122
.
Kempter
,
R.
,
Gerstner
,
W.
,
van Hemmen
,
J. L.
, &
Wagner
,
H.
(
1998
).
Extracting oscillations: Neuronal coincidence detection with noisy periodic spike input
.
Neural Comput.
,
10
,
1987
2017
.
Köppl
,
C.
(
1997
).
Phase locking to high frequencies in the auditory nerve and cochlear nucleus magnocellularis of the barn owl, Tyto alba
.
J. Neurosci.
,
17
(
9
),
3312
3321
.
Lansky
,
P.
, &
Sato
,
S.
(
1999
).
The stochastic diffusion models of nerve membrane depolarization and interspike interval generation
.
J. Peripher. Nerv. Syst.
,
4
,
27
42
.
Lapicque
,
L.
(
1907
).
Recherches quantitatives sur l'excitation électrique des nerfs traitée comme une polarisation
.
J. Physiol. (Paris)
,
9
,
620
635
.
Lindner
,
B.
,
Chacron
,
M. J.
, &
Longtin
,
A.
(
2005
).
Integrate-and-fire neurons with threshold noise: A tractable model of how interspike interval correlations affect neuronal signal transmission
.
Phys. Rev. E.
,
72
,
021911–1–21
.
Mannella
,
R.
(
2000
).
A gentle introduction to the integration of stochastic differential equations
. In
J. A. Freund & T. Pöschel
(Eds.),
Stochastic processes in physics, chemistry, and biology
(pp.
353
364
).
Berlin
:
Springer
.
Naundorff
,
B.
, &
Freund
,
J. A.
(
2002
).
Signal detection by means of phase coherence induced through phase resetting
.
Phys. Rev. E.
,
66
,
040901(R)
1
4
.
Neiman
,
A.
,
Pei
,
X.
,
Russell
,
D.
,
Wojtenek
,
W.
,
Wilkens
,
L.
,
Moss
,
F.
, et al
(
1999
).
Synchronization of the noisy electrosensitive cells in the paddlefish
.
Phys. Rev. Lett.
,
82
,
660
663
.
Ogawa
,
T.
,
Bishop
,
P. O.
, &
Levick
,
W. R.
(
1966
).
Temporal characteristics of responses to photic stimulation by single ganglion cells in the unopened eye of the cat
.
J. Neurophysiol.
,
29
,
1
30
.
Pikovsky
,
A.
,
Rosenblum
,
M.
, &
Kurths
,
J.
(
2001
).
Synchronization: A universal concept in nonlinear sciences
.
Cambridge
:
Cambridge University Press
.
Plesser
,
H.
, &
Geisel
,
T.
(
1999
).
Markov analysis of stochastic resonance in a periodically driven integrate-and-fire neuron
.
Phys. Rev. E.
,
59
,
7008
7017
.
Ratnam
,
R.
, &
Nelson
,
M. E.
(
2000
).
Nonrenewal statistics of electrosensory afferent spike trains: Implications for the detection of weak sensory signals
.
J. Neurosci.
,
20
,
6672
6683
.
Rieke
,
F.
,
Warland
,
D.
,
de Ruyter van Steveninck
,
R.
, &
Bialek
,
W.
(
1997
).
Spikes: Exploring the neural code
.
Cambridge MA
:
MIT Press
.
Rose
,
J. E.
,
Brugge
,
J. F.
,
Anderson
,
D. J.
, &
Hind
,
J. E.
(
1967
).
Phase-locked response to low-frequency tones in single auditory nerve fibers of the squirrel monkey
.
J. Neurophysiol.
,
30
,
769
793
.
Schwalger
,
T.
, &
Schimansky-Geier
,
L.
(
2008
).
Interspike interval statistics of a leaky integrate-and-fire neuron driven by gaussian noise with large correlation times
.
Phys. Rev. E.
,
77
,
031914–1-9
.
Stein
,
R.
(
1965
).
A theoretical analysis of neuronal variability
.
Biophys. J.
,
5
,
173
194
.
Talbot
,
W. H.
,
Darian-Smith
,
I.
,
Kornhuber
,
H. H.
, &
Mountcastle
,
V. B.
(
1968
).
The sense of flutter-vibration: Comparison of the human capacity with response patterns of mechanoreceptive afferents from the monkey hand
.
J. Neurophysiol.
,
31
,
301
334
.
Wilson
,
R. I.
, &
Mainen
,
Z. F.
(
2006
).
Early events in olfactory processing
.
Annu. Rev. Neurosci.
,
29
,
163
201
.