# Markov analysis of stochastic resonance in a periodically driven integrate-fire neuron

## I Introduction

Periodically modulated stochastic processes have been studied intensely over the last two decades under the paradigm of stochastic resonance: the transduction of signals is optimal in the presence of a particular amount of noise. First suggested to explain the periodicity of ice-ages , stochastic resonance has since been demonstrated in a wide range of experiments and the underlying mechanisms are well understood. A recent review of the field is given in .

The concept of stochastic resonance has met with particular attention in the neurosciences . The brain achieves an enormous signal processing performance in the presence of noise from a wide range of sources, ranging from stochastic membrane channel openings on a molecular level, via highly irregular firing patterns of individual neurons to distracting stimuli in perception. The improvement of signal transduction on all of these levels has now been demonstrated experimentally . Recently, the first direct evidence for the behavioral relevance of stochastic resonance has been reported , underlining the importance of stochastic resonance in neurobiology.

In short, neurons are threshold devices that receive an input $`I(t)`$ which charges the membrane of the neuron like a leaky capacitor. When the potential $`v(t)`$ across the membrane reaches a threshold $`\mathrm{\Theta }`$, a *spike* is fired: the membrane potential makes a brief but strong excursion (duration $`2`$ms, amplitude $`100`$mV). This spike is transmitted as output to other neurons. After the spike, the membrane potential is reset to a resting value $`v_0`$, some $`30`$mV below the threshold . As the shape of the spikes is stereotypical, information is only conveyed by the spike times.

This has led to the leaky integrate-fire model of neuronal dynamics . In between two spikes, the membrane potential is governed by

$$\tau _m\dot{v}(t)=v(t)+I(t)+\zeta (t).$$
(1)   
Here, $`\tau _m`$ is the time-constant of the membrane, which represents the internal time-scale of the neuron and $`\zeta (t)`$ is an as yet undefined noise process, comprising, e.g., stochastic membrane potential fluctuations and irregular input to the neuron from sources uncorrelated to $`I(t)`$. As the potential reaches the threshold, a spike is recorded and the potential is reset to $`v(t)=v_0`$ instantaneously.

For Gaussian white noise $`\zeta (t)`$ the evolution of the membrane potential $`v(t)`$ from reset potential to threshold is equivalent to an Ornstein–Uhlenbeck process with drift $`I(t)`$ and an absorbing boundary at $`v=\mathrm{\Theta }`$. The output of the neuron is modeled as a sequence of delta pulses $`f(t)=_k\delta (tt_k)`$ at the times of threshold crossings $`\{t_k\}=\{t|v(t)=\mathrm{\Theta }\}`$ (*spike train*). The spike train is a stochastic point process, specified entirely by the spike times $`\{t_k\}`$.

This biologically most interesting stochastic process has so far escaped a rigorous analysis, in spite of several partially successful attempts . For a list of open issues see Sec. V.C.4 of the review by Gammaitoni et al. . This is in marked contrast to the treatment of mathematically more accessible, but biologically less plausible models, such as bistable dynamic systems  and threshold devices *without reset* , in which stochastic resonance has been well established.

The essential difficulty arises from the reset after each spike: there is no well-defined membrane potential distribution for asymptotic times, as used in the case of reset-free threshold detectors. Instead we have to analyze each inter-spike-interval separately and then put these pieces together to obtain the spike train as a whole. To facilitate this, past work has assumed that the durations of all inter-spike-intervals ($`\tau _k=t_kt_{k1}`$) were identically and independently distributed (i.i.d.), i.e. that the spike train is a stationary renewal process . But in the presence of time-dependent input $`I(t)`$, this would require *identical input* within each inter-spike-interval (ISI). This is the much criticized reset assumption: If the model neuron were to describe a neuron in the auditory nerve while you are listening to a music tape, the reset assumption requires that upon the firing of each spike the tape should be rewound to exactly the position it had at the time of the last spike!

In this work we show how to analyze the response of the leaky integrate-fire neuron to periodic stimuli without undue assumptions. The distribution of the length of individual inter-spike intervals is computed numerically , and spike trains are then assembled as Markov chains from these intervals. We obtain probability distributions for the length of inter-spike intervals and the stimulus phases at which spikes occur. These distributions should be directly comparable to experiments employing sustained stimulation with periodic signals. The signal processing performance of the neuron is judged by the signal-to-noise ratio (SNR) of the output spike train. The SNR is maximal at an optimal noise amplitude for fixed stimulus frequency and at a resonance frequency for fixed noise amplitude. The latter resonance is a consequence of a time-scale matching between stimulus and membrane time-constant. All computations are verified by simulations.

In Sec. II, we show how to exploit the Markov property of the integrate-fire neuron to determine its response to sinusoidal input $`I(t)`$. The performance of the model neuron as a signal processing device is investigated in Sec. III. The results are discussed in Sec. IV.

## II Markov analysis

For an input current consisting of a constant offset and a sinusoidal component, and Gaussian white noise the Langevin equation (1) reads

$$\dot{v}(t)=v(t)+\mu +q\mathrm{cos}(\mathrm{\Omega }t+\varphi _0)+\sqrt{D}\xi (t),$$
(2)   
where time and potential have been scaled to their respective natural units $`\tau _m`$ and $`\mathrm{\Theta }`$; the reset potential is set to $`v_0=0`$. The input is characterized by the DC offset $`\mu `$, stimulus amplitude $`q`$, frequency $`\mathrm{\Omega }`$ and initial phase $`\varphi _0`$. The noise has amplitude $`\sqrt{D}`$ and autocorrelation $`\xi (t)\xi (t^{})=\delta (tt^{})`$. In the remainder of this article, we will investigate this model. For a derivation of the type of input current used here from more elementary models, see .

In the absence of noise ($`D=0`$), spikes will only be generated if

$`v_{\mathrm{}}=\underset{t\mathrm{}}{lim}v(t)=\mu +{\displaystyle \frac{q}{\sqrt{1+\mathrm{\Omega }^2}}}>1.`$

Therefore, we classify stimuli as sub-threshold if $`v_{\mathrm{}}1`$ and as supra-threshold otherwise. In this work, we will focus on the biologically more interesting sub-threshold regime . The methods presented here are applicable independent of the choice of stimulus parameters. We only require the presence of noise, i.e. $`D>0`$.

Suppose that an initial spike has occured at time $`t_0=0`$, corresponding to stimulus phase $`\varphi _0`$. The next spike follows at time $`t_1=inf\{t>t_0|v(t)1\}`$ and stimulus phase $`\varphi _1=(\mathrm{\Omega }t_1+\varphi _0)mod2\pi `$, whence $`\mathrm{cos}(\mathrm{\Omega }(tt_1)+\varphi _1)=\mathrm{cos}(\mathrm{\Omega }t+\varphi _0)`$ for $`t>t_1`$. This suggests to re-write Eq. (2) in terms of the time $`t^{}`$ that has passed since the most recent spike at phase $`\varphi `$. Thus, given this phase, the potential evolves from $`v(t^{}=0|\varphi )=v_0=0`$ until the next spike according to

$$\dot{v}(t^{}|\varphi )=v(t^{}|\varphi )+\mu +q\mathrm{cos}(\mathrm{\Omega }t^{}+\varphi )+\sqrt{D}\xi (t^{}).$$
(3)   
The next spike is fired after an interval $`\tau `$, as soon as the threshold condition is met

$$\tau =inf\{t^{}>0|v(t^{}|\varphi )1\}.$$
(4)   
The inter-spike intervals are connected by the iteration equations

$$\varphi _k=(\mathrm{\Omega }\tau _k+\varphi _{k1})mod2\pi ,t_k=t_{k1}+\tau _k,$$
(5)   
leading to the output spike train

$$f(t)=\underset{j=0}{\overset{\mathrm{}}{}}\delta (tt_j)=\underset{j=0}{\overset{\mathrm{}}{}}\delta \left(t\underset{k=1}{\overset{j}{}}\tau _k\right).$$
(6)   
The reset of the membrane potential to $`v_0=0`$ after each spike completely erases the memory of the neuron. The subsequent behavior of the neuron therefore depends on its past only through the absolute time of the spike $`t_k`$, i.e. the spike train is a Markov process.

We have thus split the task of solving the dynamics of the integrate-fire neuron into two parts. We will first solve the first-passage-time problem posed by Eqs. (3) and (4) for a given phase $`\varphi `$ of the last spike, before assembling the spike train from the inter-spike intervals according to Eqs. (5) and (6).

### A Conditional ISI distribution

The first-passage-time problem for the membrane potential posed by Eqs. (3, 4) yields the distribution $`\rho (\tau |\varphi )`$ of the inter-spike-interval lengths $`\tau `$ for a given stimulus phase $`\varphi `$ at the beginning of the interval (conditional ISI distribution). To the best of our knowledge, no analytic solution is known for this seemingly simple first-passage time problem of the Ornstein–Uhlenbeck process. The approximations suggested in  are valid in a restricted parameter range only—low stimulus frequencies in particular—and appear to yield qualitative rather than quantitative agreement with simulations.

We employ here a numerical method to compute the inter-spike-interval distributions. The method is discussed in detail in , and we only sketch it here. In the absence of an absorbing threshold the probability $`𝒫(w,t^{}|u,s^{};\varphi )`$ that the membrane potential is $`w`$ at time $`t^{}`$ if it was $`u`$ at time $`s^{}<t^{}`$ is a Gaussian distribution. The mean is given by the solution at time $`t^{}`$ of Eq. (3) for the noise-free case ($`D=0`$) with initial condition $`v(s^{}|\varphi )=u`$, while the variance is $`\sigma ^2(t^{})=\frac{D}{2}(1e^{2(t^{}s^{})})`$. Then, the inter-spike-interval distribution is given by the integral equation 

$$𝒫(1,t^{}|0,0)=_0^t^{}𝒫(1,t^{}|1,\tau )\rho (\tau |\varphi )𝑑\tau .$$
(7)   
This equation is solved for $`\rho `$ using standard techniques . Source code is available on request.

As shown in Fig. 1, the conditional inter-spike-interval distributions $`\rho (\tau |\varphi )`$ may depend strongly on $`\varphi `$. First, they contain a series of exponentially decaying peaks that are separated by the stimulus period $`T=2\pi /\mathrm{\Omega }`$. These peaks represent spikes that are well phase-locked to the stimulus and we will refer to them as *periodic peaks*. An additional peak appears at short intervals $`\tau `$ for certain phases $`\varphi `$. This peak reflects the rise time of the membrane potential towards threshold. Its location is not related to the stimulus period $`T`$, but reflects the intrinsic time-scale of the neuron and defines its refractory time, i.e. the minimum interval between two spikes. Thus, we will refer to this peak as the *refractory peak*. It corresponds to two or more spikes fired in rapid succession within a single stimulus period (a *burst*). There is thus a qualitative dependence of the distributions $`\rho `$ on the phase $`\varphi `$ that can lead to interesting consequences for the firing behavior of the neuron.

The relation between periodic and refractory peaks depends on the stimulus parameters, particularly on the frequency and the noise amplitude. We will discuss this relationship in Sec. II C.

### B Markov process in phase

Let us now turn to the problem of assembling spike trains from inter-spike-intervals according to Eqs. (5, 6). The length of an interval following a spike at time $`t`$ and stimulus phase $`\varphi =[\mathrm{\Omega }t+\varphi _0]mod2\pi `$ is distributed according to $`\rho (\tau |\varphi )`$. Therefore, the probability that the next spike will occur at phase $`\psi `$ is given by

$$𝒯\left(\psi |\varphi \right)=_0^{\mathrm{}}\rho (\tau |\varphi )\delta (\psi [\mathrm{\Omega }\tau +\varphi ]mod2\pi )\frac{d\tau }{\mathrm{\Omega }}.$$
(8)   
We will call $`𝒯\left(\psi |\varphi \right)`$ the transition probability of the spike phase. We will now consider the Markov process of the spike phases $`\varphi _k`$ instead of the Markov process made up of the spike times $`t_k`$.

If we define the spike phase distribution $`\chi ^{(k)}\left(\varphi \right)`$ as the probability (across an ensemble of neurons or repetitions of an experiment) that the $`k^{\text{th}}`$ spike in a train will be fired at stimulus phase $`\varphi `$, then this probability will evolve according to

$$\chi ^{(k+1)}\left(\psi \right)=_0^{2\pi }𝒯\left(\psi |\varphi \right)\chi ^{(k)}\left(\varphi \right)𝑑\varphi .$$
(9)   
As the neuron fires repetitively while driven by a stationary periodic stimulus, the spike train emitted by the neuron will approach a stationary Markov process with phase distribution

$$\chi ^{(s)}\left(\psi \right)=\underset{k\mathrm{}}{lim}\chi ^{(k)}\left(\psi \right)=_0^{2\pi }𝒯\left(\psi |\varphi \right)\chi ^{(s)}\left(\psi \right)𝑑\varphi .$$
(10)   
The stationary phase distribution $`\chi ^{(s)}\left(\psi \right)`$ is the eigenfunction to eigenvalue $`1`$ of the kernel $`𝒯\left(\psi |\varphi \right)`$, and is guaranteed to exist because this kernel is a conditional probability distribution . Any initial phase distribution will converge to the unique stationary solution provided that $`𝒯\left(\psi |\varphi \right)>0`$ everywhere . That the latter condition holds in the presence of noise can be seen as follows. For sub-threshold stimuli, noise may drive the potential across the firing threshold at any time $`\tau >0`$ in principle, yielding a possibly tiny, but non-zero probability of spikes at any phase. The same argument holds true for supra-threshold stimuli, where noise may keep potential below threshold up to any time. In the absence of noise, neither convergence nor uniqueness are assured.

To facilitate numerical treatment, we discretize the phase. Since the conditional inter-spike-interval distributions $`\rho (\tau |\varphi )`$ are smooth in both time and phase due to the presence of noise in the input, this discretization will introduce only minor numerical errors. It is largely equivalent to applying numerical methods to solve the kernel eigenvalue problem . Using $`L`$ bins of width $`\mathrm{\Delta }\psi `$ ($`\mathrm{\Delta }\psi =2\pi /L`$) we obtain the spike phase distribution vector

$$𝝌=(\chi _0,\chi _1,\mathrm{},\chi _{L1})^{\text{tr}},\chi _j=_{j\mathrm{\Delta }\psi }^{(j+1)\mathrm{\Delta }\psi }\chi (\psi )𝑑\psi ,$$
(11)   
and the phase transition matrix $`𝐓`$ with elements

$$𝐓_{jk}=_{j\mathrm{\Delta }\psi }^{(j+1)\mathrm{\Delta }\psi }𝒯\left(\psi |k\mathrm{\Delta }\psi \right)𝑑\psi ,j,k=0,\mathrm{},L1.$$
(12)   
The evolution equation (9) simplifies from convolution to matrix-vector multiplication

$$𝝌^{(k+1)}=𝐓𝝌^{(k)},$$
(13)   
and the stationary distribution $`𝝌^{(s)}`$ is the eigenvector to eigenvalue $`1`$ of the matrix $`𝐓`$. We have thus reduced the Markov process to a Markov chain.

In practice, we obtain the transition matrix $`𝐓`$ by numerically evaluating equations (8) and (12), with $`\rho (\tau |\varphi )`$ from Eq. (7). The stationary distribution is then found using standard eigenvector routines. For all data shown here, we used the discretization $`L=72`$, $`\mathrm{\Delta }\psi =\pi /36=5^{}`$. In figures of transition matrices and phase distributions the axis will run from $`\pi `$ to $`\pi `$ as this renders structures more clearly.

An example for the phase evolution of an initially uniform distribution towards the stationary state under the influence of a transition matrix $`𝐓`$ is given in Fig. 2. To “read” the transition matrix, note that the matrix columns correspond to the phase $`\varphi _k`$ of the spike preceding the interval, the rows to the phase $`\varphi _{k+1}`$ of the spike terminating it. The phase axes run from $`\pi `$ to $`\pi `$ from bottom to top in phase distribution vectors $`𝝌`$ and the rows of the transition matrix $`𝐓`$, and *from right to left* across the columns of $`𝐓`$. Thus, the horizontal bar in the transition matrix shown in Fig. 2 indicates that for most values of $`\varphi _k`$ the next spike will occur around $`\varphi _{k+1}\pi /6`$. This bar corresponds to the periodic peaks of the ISI distributions. For $`\pi /4\varphi _k\pi /6`$, the matrix is dominated by a “finger”, running parallel to the matrix diagonal. Within this range of phases, a spike will be followed by another spike at a slightly later phase, as shown in Fig. 2b. Figuratively speaking, the neuron fires a burst of spikes, but there is always a chance that two subsequent spikes will be one or more stimulus periods apart, even though they are close in phase: in the Markov chain description, all information about actual interval lengths is lost. The finger results from the refractory peak of the ISI distributions.

Figure 3 shows the dependence of transition matrix and stationary phase distribution on the noise amplitude for slow stimuli ($`T10`$). For low noise, the transition matrix is dominated by the horizontal bar, which intersects with the matrix diagonal, indicating a stochastic fixed point. This results in a sharply peaked spike phase distribution. At intermediate noise, the finger is more pronounced, while the bar barely touches the matrix diagonal, leading to a stochastic limit cycle with two preferred phases: the neuron often fires bursts of two successive spikes. At high noise, the finger stretches all along the matrix diagonal, while the horizontal bar has disappeared altogether. The neuron fires rapidly, but largely uncorrelated with the stimulus and the phase distribution is virtually flat.

This means that for very low noise the spike train of the neuron is nearly a stationary renewal process with inter-spike-intervals i.i.d. according to $`\rho (\tau |\psi ^{})`$. Here $`\psi ^{}`$ is the location of the maximum of the stationary phase distribution, which depends not only on the stimulus parameters, but also on the noise amplitude. For high noise amplitudes, the response of the neuron is largely independent of the stimulus, and may thus be described by a stationary renewal process as well—the ISIs reduce to the refractory peak. But at intermediate noise levels—i.e. those essential to the observation of stochastic resonance—the stationary phase distribution may be multimodal. Thus the correlations between the phases of subsequent spikes have to be taken into account using the Markov ansatz. Multimodal phase distributions as discussed here are not just hypothetical: they have been observed in sensory neurons of goldfish upon stimulation with sinusoidal water waves .

For fast stimuli ($`T10`$), the stationary phase distribution smears out much more along the phase axis, and does not show multimodality, because the refractory time of the neuron becomes comparable to the stimulus period and bursting is no longer possible, see Fig. 4. At low to intermediate noise, the distribution is too wide to be replaced by its mode as in the renewal ansatz, but still sufficiently narrow to provide for a response that is well phase-locked to the stimulus. Therefore, the Markov approach is essential for high frequency stimuli as well.

### C Stationary ISI distribution

Once the stationary phase distribution is known, the inter-spike-interval distribution of the stationary firing process is obtained by averaging the conditional ISI distributions over phase

$$\rho (\tau )=_0^{2\pi }\rho (\tau |\psi )\chi ^{(s)}\left(\psi \right)𝑑\psi .$$
(14)   
The average interval length thus is

$$\tau =_0^{\mathrm{}}\tau \rho (\tau )𝑑\tau .$$
(15)   
$`\rho (\tau )`$ is the inter-spike-interval distribution that we expect to find in experiments with tonic stimulation. In contrast to a stationary renewal process, this averaged ISI distribution does *not* contain a full description of the spike train.

Typical ISI distributions $`\rho (\tau )`$ are given in Figs. 5 and 6 for the same parameters as used in Figs. 34, respectively. For low noise, they contain only periodic peaks, located precisely at integer multiples of the stimulus period $`T`$: the neuron can only fire in a small time window within each period, and several periods may be skipped in between spikes. This indicates a firing pattern that is well phase-locked to the stimulus. ISI distributions with comparable structure have been found in neurons of the auditory system in different species . For high noise, the ISI distributions reduce to the refractory peak, i.e. a largely random firing pattern.

For intermediate noise, the ISI distributions depend strongly on the stimulus frequency. For high frequency (Fig. 6), we find merely a superposition of periodic and refractory peaks: spikes preferentially occur at intervals that are multiples of the stimulus period, but this phase-locking is weak. This is very different for slow stimuli (Fig. 5), where the refractory peak is clearly separated from a wide peak at $`\tau =T=40`$, the latter exposing some sub-structure. This can be understood as follows. The maximum of $`\rho (\tau )`$ at $`\tau =T`$ corresponds to two spikes fired each at the optimal phase in two subsequent periods. In contrast, if a period that contained a burst of two spikes is followed by another period containing a burst, then typically the first spike will be slightly earlier than the optimal phase, the second one a bit later. Thus, the interval between the second spike of the first burst and the first spike of the second burst is shorter than the stimulus period, leading to the side-peak at $`\tau 35`$. The bursts themselves give rise to the refractory peak. This again indicates that the spike train is not a stationary renewal process.

Along with results obtained using the Markov chain approach, Figs. 36 display phase and ISI distributions obtained from simulated trains of 20,000 spikes. The agreement between Markov model and simulation is excellent. Source code for the simulation based on  is available on request.

## III Stochastic Resonance

To assess the performance of the integrate-fire neuron as a signal processing device, we evaluate the signal-to-noise ratio (SNR) of the spike train generated in response to periodic input. In doing so, one should keep in mind the purpose of the output spike train. It has to convey information to other neurons in the brain *within a certain time window*, as the brain has to respond quickly to stimuli. Therefore, the relevant quantity is the signal-to-noise ratio that can be achieved by measuring the spike train over a finite observation time $`T_o`$ .

### A Signal-to-noise ratio

The one-sided power spectral density of a stationary spike train $`f(t)`$ \[as defined in Eq. (6)\] over a time interval $`T_o`$ is 

$`S_{T_o}^{}(\omega )`$ $`={\displaystyle \frac{1}{\pi T_o}}\left|{\displaystyle _0^{T_o}}f(t)e^{i\omega t}𝑑t\right|^2`$ (17)   
$`={\displaystyle \frac{1}{\pi T_o}}{\displaystyle \underset{j,k}{\overset{t_j,t_k<T_o}{}}}e^{i\omega (t_jt_k)}.`$   
The average is to be taken over the ensemble of all spike trains, that is, over the set of all conditional ISI distributions and their $`(jk)`$-fold convolutions. This problem appears intractable.

The situation is greatly simplified if $`\omega `$ is the stimulus frequency $`\mathrm{\Omega }`$ or one of its harmonics. Expressing the spike times as $`t_j=(m_j+\frac{\psi _j}{2\pi })T`$, Eq. (17) for $`\omega =n\mathrm{\Omega }`$ simplifies to

$$S_{T_o}^{}(n\mathrm{\Omega })=\frac{1}{\pi T_o}\underset{j,k}{\overset{t_j,t_k<T_o}{}}e^{in(\psi _j\psi _k)},$$
(18)   
where $`n,m_j`$ are integers, $`\psi _j[\mathrm{\hspace{0.17em}0},2\pi )`$, and $`T=2\pi /\mathrm{\Omega }`$ is the stimulus period. In the observation period $`T_o`$, on average $`T_o/\tau `$ spikes will occur, regardless of the detailed structure of the spike train. We therefore fix the upper limit of the summation at $`N_o=T_o/\tau `$, where $`\tau `$ is the average interval length from Eq. (15) and $`x`$ is the largest integer not exceeding $`x`$. This yields as an approximation

$$S_{T_o}^{}(n\mathrm{\Omega })S_{T_o}(n\mathrm{\Omega })=\frac{1}{\pi N_o\tau }\underset{j,k=1}{\overset{N_o}{}}e^{in(\psi _j\psi _k)}.$$
(19)   
The task of computing an expectation with respect to all possible spike trains is now reduced to that of averaging over all possible sequences of *spike phases*. Their distribution and correlations are completely characterized by the transition matrix $`𝐓`$, permitting for evaluation of Eq. (19) in closed form. The actual calculation is straightforward albeit lengthy algebra and is provided in the appendix. The final result may be written as

$$S_{T_o}(n\mathrm{\Omega })=\frac{1}{\pi \tau }\left[\mathrm{\hspace{0.17em}1}+A(n,N_o)+(N_o1)B(n)\right]$$
(20)   
where the functions $`A(n,N_0)`$ and $`B(n)`$ are given in the appendix. Note that $`A(n,N_o)`$ is bounded as $`N_o\mathrm{}`$. For a Poissonian spike train, both $`A`$ and $`B`$ are identically zero, yielding a white power spectrum .

At first, it might seem surprising that the spectrum contains a term, $`(N_o1)B(n)`$, that scales linearly with the number of spikes in the train. This is a consequence of the periodic component of the spike train introduced by the driving stimulus, leading to a mixed spectrum consisting of a continuous background and a discrete spectrum of harmonics . For infinite observation time, i.e. $`N_o\mathrm{}`$, this gives rise to the terms $`\delta (\omega n\mathrm{\Omega })`$ in the power spectrum.

A typical power spectrum is shown in Fig. 7, indicating close agreement of Eq. (20) with results obtained by numerical Fourier transformation of simulated spike trains. The approximation made in fixing the summation limit in Eq. (19) is therefore well justified. The dip in the noise background of the spectrum at low frequencies is a consequence of the refractory period of the neuron, while the weak hump at $`\omega 1`$ indicates the presence of bursts . Spectra consisting only of this background have been found in neurons of higher cortical areas of monkeys in the absence of periodic input .

Since the power spectral density can only be evaluated in closed form at multiples of the stimulus frequency, we approximate the noise background as Poissonian white noise $`S_P=(\pi \tau )^1`$ of a spike train of equal intensity . The signal-to-noise ratio obtainable from the spike train within the observation time $`T_o`$ is therefore given by

$$\text{SNR}_{T_o}=\frac{S_{T_o}(\mathrm{\Omega })}{S_P}=1+A(1,\frac{T_o}{\tau })+(\frac{T_o}{\tau }1)B(1).$$
(21)   
The signal-to-noise ratio for three different stimulus frequencies is shown in Fig. 8 vs. the noise amplitude, again in excellent agreement with simulation results. Stochastic resonance (SR) is clearly present at all frequencies, as the SNR attains its maximum for an intermediate noise level. The striking new feature is that the overall maximum in the SNR is reached at an intermediate frequency $`\mathrm{\Omega }_r\pi /3`$, which we thus call the *resonance frequency*. The same qualitative dependence of the SNR on noise amplitude and stimulus frequencies is observed over a wide range of stimulus parameters, including weakly supra-threshold cases ($`0.4\mu <1`$, $`0.4q/(1\mu )1.2`$; data not shown).

Note that the stochastic resonance reported in an earlier paper  is an artifact of the renewal ansatz employed in that work. There, the stimulus phase is reset to an arbitrarily chosen value $`\varphi _0`$ after each spike, and the signal-to-noise ratio is computed for an infinite observation time. The SNR is maximized for that noise level at which the periodic peaks of the ISI distribution $`\rho (\tau |\varphi _0)`$ are centered about the multiples of the stimulus period $`T`$. But if, for low noise, one uses for each noise level $`D`$ a different $`\varphi _0(D)`$, namely the mode of the stationary phase distribution as discussed in Sec. II B, the periodic peaks are at multiples of $`T`$ for all noise intensities, whence the SNR does not drop off for $`D0`$ and no resonance occurs (data not shown). This observation underlines the importance of the Markov approach.

### B Time-scale matching

In contrast to stochastic resonance in dynamical systems, SR with respect to the noise amplitude is not induced by the matching of time-scales in threshold systems, but results from stochastic linearization of the response function of the neuron . In contrast, the additional resonance along the frequency axis arises in the integrate-fire neuron as a consequence of matching the stimulus period to the intrinsic time scale of the neuron in an appropriate manner. This is demonstrated in Fig. 9. For a stimulus at the resonance frequency $`\mathrm{\Omega }_r`$, the peak at $`\tau =T`$ in the stationary ISI distribution can “grow” in place as noise is increased, without being disturbed by the refractory peak. Indeed, the latter arises at the location of the first periodic peak and shifts away from $`\tau =T`$ only for very large noise. In this way, the firing rate of the neuron can be increased without loosing the phase-locking to the stimulus. Compare this to the cases of lower (Fig. 5) and higher (Fig. 6) frequencies: in both cases, high firing rates can only be achieved by raising the noise amplitude to a point where the refractory peak has either replaced ($`\mathrm{\Omega }<\mathrm{\Omega }_r`$) or smeared out ($`\mathrm{\Omega }>\mathrm{\Omega }_r`$) the periodic peaks, resulting in a firing pattern poorly phase-locked to the stimulus.

This competition of precision and intensity is demonstrated by a phenomenological ansatz for the SNR. A measure of phase-locking between stimulus and response is the vector strength $`C_s=\left|e^{i\psi }\right|`$, where $`\psi `$ are the spike phases . $`C_s=1`$ indicates perfect and $`C_s=0`$ no locking. If the neuron attempts to measure the degree of phase-locking from a train of $`N=T_o/\tau `$ spikes, the quality of measurement will be $`\sqrt{N}`$. Thus, we expect that the signal-to-noise ratio will roughly given by

$$\text{SNR}_{\text{phen}}C_s\sqrt{N}=C_s\sqrt{\frac{T_o}{\tau }}.$$
(22)   
Figure 10 demonstrates that this simple model describes the behavior of the SNR well. In particular, the two-fold stochastic resonance is reproduced.

In short, to elicit a strong output signal from the model neuron, a sufficient input noise level is required. But this comes at a cost, as the quality of the output, i.e. the precision of the phase locking, deteriorates as noise is added. The maximum SNR represents the optimal compromise between signal strength and quality.

## IV Discussion

In this paper, we have shown that the periodically driven integrate-fire neuron can be analyzed in the framework of a Markov process. This avoids the unrealistic assumption of a stimulus reset after each spike, the most serious shortcoming of previous work , and this answers question (1) raised by Gammaitoni et al. in Sec. V.4.C of their review . Their second questions concerns the fact that the neural membrane is a rectifier: even a strong negative input current will not lower the membrane potential more than a few millivolts below the reset potential $`v_0`$. This would indeed be a problem if the DC offset $`\mu `$ of the input were much smaller than the amplitude $`q`$ of the AC stimulus. Preliminary evidence suggests that the best fit of inter-spike-interval distributions generated by the model with experimental data from the cat’s auditory system  is obtained for sub-threshold stimuli with $`\mu q`$. In this regime, the membrane potential is quickly raised to $`v_0+\mu `$ and then oscillates around this level, unaffected by rectification. Finally, Gammaitoni and co-authors question the validity of the approximations used to compute the ISI distributions in . This matter is avoided here by numerically computing these distributions. A study of the validity of approximate closed-form ISI distributions will be given elsewhere .

The Markov formalism presented in this paper is applicable to any periodically driven stochastic process with a reset. The only required ingredients are the conditional first-passage-time distributions $`\rho (\tau |\varphi )`$ and the iteration equations (5). The generalization to more complex stimuli, e.g. including amplitude modulation, is straightforward.

With the Markov machinery at hand, we have demonstrated that the signal-to-noise ratio of the output of the neuron is maximized at an optimal noise amplitude for fixed frequency and at a resonance frequency for fixed noise intensity. Stochastic resonance with respect to the stimulus frequency, termed *bona fide* stochastic resonance, has been described in bistable systems before . Therefore, our findings for a non-dynamical threshold neuron extend the universality of stochastic resonance to the case of *bona fide* SR. Recent criticism  of the original definition of *bona fide* SR, based on residence time distributions, does not apply to our study.

Neurons in the auditory system can phase lock to acoustic stimuli with high acuity and utilize this for the precise localization of sound sources . Our results show that strong signals that are well phase locked to a stimulus may be achieved in spite of the noise ubiquitous in the neural system. Stochastic resonance might therefore be one of the underlying mechanisms of stereo hearing. First qualitative comparisons indicate good agreement between response properties of the integrate-fire neuron and of auditory neurons. An intriguing question in this respect is the relevance of the *bona fide* SR to the neural system. It may serve to tune neurons as bandpass filters of a special kind: only stimuli in a certain frequency window will be transmitted with high intensity *and* precise phase locking. A detailed study will be the topic of a future publication.

###### Acknowledgements.

This work was supported by Deutsche Forschungsgemeinschaft through SFB 185 “Nichtlineare Dynamik”. HEP gratefully acknowledges the hospitality of the Laboratory for Neural Modeling, Frontier Research Program, RIKEN, Wako-shi, Saitama, Japan, where this work started.    

## A Computing the power spectral density

To prove Eq. (20), i.e.

$`S_{T_o}(n\mathrm{\Omega })`$ $`={\displaystyle \frac{1}{\pi M\tau }}{\displaystyle \underset{j,k=1}{\overset{M}{}}}e^{in(\psi _j\psi _k)}`$   
$`={\displaystyle \frac{1}{\pi \tau }}\left[\mathrm{\hspace{0.17em}1}+A(n,M)+(M1)B(n)\right],`$   
we split the double sum into the diagonal and off-diagonal terms

$$S_{T_o}(n\mathrm{\Omega })=\frac{1}{\pi \tau }\left[1+h_M(n\mathrm{\Omega })+h_M^{}(n\mathrm{\Omega })\right],$$
(A1)   
$$h_M(n\mathrm{\Omega })=\frac{1}{M}\underset{k=1}{\overset{M}{}}\underset{j=1}{\overset{Mk}{}}e^{in(\psi _{k+j}\psi _k)},$$
(A2)   
the asterisk denoting complex conjugation and $`M=T_o/\tau `$.

Since we are considering a stationary Markov process, all $`\psi _k`$ are identically distributed according to $`𝝌^{(s)}`$, while correlations between $`\psi _k`$ and $`\psi _{k+j}`$ are given by the $`j^{\mathrm{th}}`$ power of the transition matrix $`𝐓`$ yielding

$$e^{in(\psi _{k+j}\psi _k)}=\widehat{𝐚}(n)^{\mathrm{tr}}𝐓^j\widehat{𝐛}(n)$$
(A3)   
with vectors

$`\widehat{𝐚}^{\mathrm{tr}}(n)=(1,e^{in\mathrm{\Delta }\psi },e^{2in\mathrm{\Delta }\psi },\mathrm{},e^{(L1)in\mathrm{\Delta }\psi })`$

$`\widehat{𝐛}^{\mathrm{tr}}(n)=(𝝌^{(s)}(0),\mathrm{},e^{(L1)in\mathrm{\Delta }\psi }𝝌^{(s)}((L1)\mathrm{\Delta }\psi ))`$

Upon inserting Eq. (A3) into Eq. (A2), we observe that the expression for $`h_M`$ depends only on $`j`$ but not on $`k`$ so that we may perform the outer summation to obtain

$`h_M(n\mathrm{\Omega })=\widehat{𝐚}^{\mathrm{tr}}(n)\left[{\displaystyle \frac{1}{M}}{\displaystyle \underset{j=1}{\overset{M1}{}}}(Mj)𝐓^j\right]\widehat{𝐛}(n).`$

Diagonalizing $`𝐓`$ leads to

$$h_M(n\mathrm{\Omega })=𝐚(n)^{\text{tr}}𝐒^{(M)}𝐛(n).$$
(A4)   
Here, the diagonal matrix $`𝐒^{(M)}`$ is given by

$$S_{mm}^{(M)}=\{\begin{array}{cc}\frac{M1}{2}\hfill & \text{for }m=1,\hfill \\ \frac{\lambda _m}{1\lambda _m}+\frac{1}{M}\frac{\lambda _m(\lambda _m^M1)}{(\lambda _m^M1)^2}\hfill & \text{for}m>1,\hfill \end{array}$$
(A5)   
with

$`𝐓=𝐂𝐋𝐂^1,𝐋=\mathrm{diag}(1>|\lambda _2|\mathrm{}|\lambda _L|),`$

$`𝐚(n)=𝐂^{\text{tr}}\widehat{𝐚}(n),𝐛(n)=𝐂^1\widehat{𝐛}(n).`$

Inserting Eq. (A4) into Eq. (A1), we have

$$S_{T_o}(n\mathrm{\Omega })=\frac{1}{\pi \tau }\left[1+2\text{Re}\left(𝐚^{\text{tr}}(n)𝐒^{(M)}𝐛(n)\right)\right].$$
(A6)   
Finally, we split the matrix $`𝐒^{(M)}`$ into the parts pertaining to the discrete and the continuous parts of the spectrum and define the functions $`A`$ and $`B`$

$`𝐒^{(M)}=\text{diag}(\frac{M1}{2},0,\mathrm{},0)+\text{diag}(0,S_{22}^{(M)},\mathrm{},S_{LL}^{(M)}),`$

$`A(n,M)=2\text{Re}\left[𝐚^{\text{tr}}(n)\text{diag}(0,S_{22}^{(M)},\mathrm{},S_{LL}^{(M)})𝐛(n)\right],`$

$`B(n)=\text{Re}\left[a_1(n)b_1(n)\right].`$

Rewriting Eq. (A6) accordingly, we arrive at the desired expression for the power spectral density

$`S_{T_o}(n\mathrm{\Omega })={\displaystyle \frac{1}{\pi \tau }}\left[\mathrm{\hspace{0.17em}1}+A(n,M)+(M1)B(n)\right].`$

To see that $`A`$ is bounded as $`M\mathrm{}`$, note that $`A`$ depends on $`M`$ only through the diagonal entries of $`𝐒^{(M)}`$ with $`m>1`$ and $`|\lambda _m|<1`$. For these we have

$`\underset{M\mathrm{}}{lim}|S_{mm}^{(M)}|=\left|{\displaystyle \frac{\lambda _m}{1\lambda _m}}\right|<\mathrm{},m>1.`$