Next: 5.3 Escape noise Up: 5. Noise in Spiking Previous: 5.1 Spike train variability

Subsections

# 5.2 Statistics of spike trains

In this section, we introduce some important concepts for the statistical description of neuronal spike trains. A central notion will be the interspike interval distribution which is discussed in the framework of a generalized input-dependent renewal theory. We start in Section 5.2.1 with the definition of renewal systems and turn then in Section 5.2.2 to interval distributions. The relation between interval distributions and neuronal models will be the topic of Sections 5.3 and 5.5.

## 5.2.1 Input-dependent renewal systems

We consider a single neuron such as an integrate-and-fire or SRM unit. Let us suppose that we know the last firing time < t of the neuron and its input current I. In formal spiking neuron models such as the SRM, the membrane potential u is then completely determined, i.e.,

 u(t|) = (t - ) + (t - , s) I(t - s) ds , (5.1)

cf. Eq. (4.24). In particular, for the integrate-and-fire model with membrane time constant and capacity C we have

 u(t|) = ur exp -   +  exp I(t - s) ds , (5.2)

cf. Eq. (4.10). In general, part or all of the input current I could arise from presynaptic spikes. Here we simply assume that the input current I is a known function of time.

Given the input and the firing time we would like to predict the next action potential. In the absence of noise, the next firing time t(f) of a neuron with membrane potential (5.1) is determined by the threshold condition u = . The first threshold crossing occurs at

 t(f) = mint >  | u(t|) . (5.3)

In the presence of noise, however, we are no longer able to predict the exact firing time of the next spike, but only the probability that a spike occurs. The calculation of the probability distribution of the next firing time for arbitrary time-dependent input I is one of the major goals in the theory of noisy spiking neurons.

Equations. (5.1) and (5.2) combined with a (stochastic) spike generation procedure are examples of input-dependent renewal systems. Renewal processes are a class of stochastic point processes that describe a sequence of events (spikes) in time (Cox, 1962; Papoulis, 1991). Renewal systems in the narrow sense (stationary renewal processes), presuppose stationary input and are defined by the fact that the state of the system, and hence the probability of generating the next event, depends only on the age' t - of the system, i.e., the time that has passed since the last event (last spike). The central assumption of renewal theory is that the state does not depend on earlier events (i.e., earlier spikes of the same neuron). The aim of renewal theory is to predict the probability of the next event given the age of the system.

Here we use the renewal concept in a broader sense and define a renewal process as a system where the state at time t, (and hence the probability of generating an event at t), depends both on the time that has passed since the last event (i.e., the firing time ) and the input I(t'), < t' < t, that the system received since the last event. Input-dependent renewal systems are also called modulated renewal processes (Reich et al., 1998), non-stationary renewal systems (Gerstner, 1995,2000b), or inhomogeneous Markov interval processes (Kass and Ventura, 2001). The aim of a theory of input-dependent renewal systems is to predict the probability of the next event, given the timing of the last event and the input I(t') for < t' < t.

### 5.2.1.1 Example: Light bulb failure as a renewal system

A generic example of a (potentially input-dependent) renewal system is a light bulb. The event is the failure of the bulb and its subsequent exchange. Obviously, the state of the system only depends on the age of the current bulb, and not on that of any previous bulb that has already been exchanged. If the usage pattern of the bulbs is stationary (e.g., the bulb is switched on during 10 hours each night) then we have a stationary renewal process. If usage is irregular (higher usage in winter than in summer, no usage during vacation), the aging of the bulb will be more rapid or slower depending on how often it is switched on and off. We can use input-dependent renewal theory if we keep track of all the times we have turned the switch. The input in this case are the switching times. The aim of renewal theory is to calculate the probability of the next failure given the age of the bulb and the switching pattern.

## 5.2.2 Interval distribution

The estimation of interspike interval (ISI) distributions from experimental data is a common method to study neuronal variability given a certain stationary input. In a typical experiment, the spike train of a single neuron (e.g., a neuron in visual cortex) is recorded while driven by a constant stimulus. The stimulus might be an external input applied to the system (e.g., a visual contrast grating moving at constant speed); or it may be an intracellularly applied constant driving current. The spike train is analyzed and the distribution of intervals sk between two subsequent spikes is plotted in a histogram. For a sufficiently long spike train, the histogram provides a good estimate of the ISI distribution which we denote as P0(s); cf. Fig. 5.1A. We will return to the special case of stationary input in subsection 5.2.4.

We now generalize the concept of interval distributions to time-dependent input. We concentrate on a single neuron which is stimulated by a known input current I(t) and some unknown noise source. We suppose that the last spike occurred at time and ask the following question. What is the probability that the next spike occurs between t and t + t, given the spike at and the input I(t') for t' < t? For t 0, the answer is given by the probability density of firing PI(t|). Hence, PI(t|) dt is the probability to find a spike in the segment [t1, t2], given that the last spike was at < t1. The normalization of PI(t|) is

 PI(t | ) dt = 1 - pinactI (5.4)

where pinactI denotes the probability that the neurons stays inactive and will never fire again. For excitatory input and a sufficient amount of noise the neuron will always emit further spikes at some point. We therefore assume in the following that pinactI vanishes.

The lower index I of PI(t|) is intended to remind us that the probability density PI(t|) depends on the time course of the input I(t') for t' < t. Since PI(t|) is conditioned on the spike at , it can be called a spike-triggered spike density. We interpret PI(t | ) as the distribution of interspike intervals in the presence of an input current I. In the following, we will refer to PI as the input-dependent interval distribution; see Fig. 5.1B. For renewal systems with stationary input PI(t|) reduces to P0(t - ).

## 5.2.3 Survivor function and hazard

The interval distribution PI(t|) as defined above is a probability density. Thus, integration of PI(t|) over time yields a probability. For example, PI(t'|) dt' is the probability that a neuron which has emitted a spike at fires the next action potential between and t. Thus

 SI(t|) = 1 - PI(t'|) dt' (5.5)

is the probability that the neuron stays quiescent between and t. SI(t|) is called the survivor function: it gives the probability that the neuron survives' from to t without firing.

The survivor function SI(t|) has an initial value SI(|) = 1 and decreases to zero for t. The rate of decay of SI(t|) will be denoted by (t|) and is defined by

 (t|) = - SI(t|) / SI(t|) . (5.6)

In the language of renewal theory, (t|) is called the age-dependent death rate' or hazard' (Cox, 1962; Cox and Lewis, 1966).

Integration of Eq. (5.6) yields the survivor function

 SI(t|) = exp - (t'|) dt' . (5.7)

According to the definition of the survivor function in Eq. (5.5), the interval distribution is given by

 PI(t|) = - SI(t|) = (t|) SI(t|) , (5.8)

which has a nice intuitive interpretation: In order to emit its next spike at t, the neuron has to survive the interval (, t) without firing and then fire at t. The survival probability is SI(t|) and the hazard of firing a spike at time t is (t|) which explains the two factors on the right-hand side of Eq. (5.8). Inserting Eq. (5.7) in (5.8), we obtain an explicit expression for the interval distribution in terms of the hazard:

 PI(t|) = (t|) exp - (t'|) dt' . (5.9)

On the other hand, given the interval distribution we can derive the hazard from

 (t|) = - = -  . (5.10)

Thus, each of the three quantities (t|), PI(t|), and SI(t|) is sufficient to describe the statistical properties of an input-dependent renewal system. For stationary renewal systems, Eqs. (5.5)-(5.10) hold with the replacement
 PI(t|) P0(t - ) (5.11) SI(t|) S0(t - ) (5.12) (t|) (t - ) . (5.13)

Eqs. (5.5) - (5.10) are standard results of renewal theory (Perkel et al., 1967a; Gerstein and Perkel, 1972; Cox, 1962; Perkel et al., 1967b; Cox and Lewis, 1966).

### 5.2.3.1 Example: From interval distribution to hazard function

Let us suppose that we have found under stationary experimental conditions an interval distribution that can be approximated as

 P0(s) = (5.14)

with a constant a0 > 0; cf. Fig. 5.2A. From Eq. (5.10), the hazard is found to be

 (s) = (5.15)

Thus, during an interval after each spike the hazard vanishes. We may interpret as the absolute refractory time of the neuron. For s > the hazard increases linearly, i.e., the longer the neuron waits the higher its probability of firing. In Section 5.3, the hazard (5.15) will be motivated by a non-leaky integrate-and-fire neuron subject to noise.

### 5.2.3.2 Example: From hazard functions to interval distributions

Interval distributions and hazard functions have been measured in many experiments. For example, in auditory neurons of the cat driven by stationary stimuli, the hazard function (t - ) increases, after an absolute refractory time, to a constant level (Goldberg et al., 1964). We approximate the time course of the hazard function as

 (s) = (5.16)

with parameters ,, and ; Fig. 5.2B. In Section 5.3 we will see how the hazard (5.16) can be related to neuronal dynamics. Given the hazard function, we can calculate the survivor function and interval distributions. Application of Eq. (5.7) yields

 S0(s) = (5.17)

The interval distribution is given by P0(s) = (sS0(s). Interval distribution, survivor function, and hazard are shown in Fig. 5.2B.

### 5.2.3.3 Example: Poisson process

Let us compare the hazard functions of the two previous examples to the hazard of a homogeneous Poisson process that generates spikes stochastically at a fixed rate . Since different spikes are independent, the hazard of a Poisson process is constant (s) . In particular, there is no dependence of the hazard upon the last or any earlier spike. From Eq. (5.8) we find the survivor function S0(s) = exp[-  s]. The interval distribution is exponential

 P0(s) =  e- s    for s > 0 . (5.18)

Interval distribution and survivor function of a Poisson neuron with constant rate are plotted in Fig. 5.3A. The most striking feature of Fig. 5.3A is that the interval distribution has its maximum at s = 0 so that extremely short intervals are most likely. In contrast to a Poisson process, real neurons show refractoriness so that the interval distribution P0(s) vanishes for s 0

A simple modification of the Poisson process allows us to incorporate absolute refractoriness. We define a hazard function

 (s) =  . (5.19)

We call a process with hazard function (5.19) a Poisson neuron with absolute refractoriness. It generates a spike train with an interval distribution

 P0(s) =  ; (5.20)

see Fig. 5.3B. We may compare the hazard function of the Poisson neuron with absolute refractoriness with the more realistic hazard of Eq. (5.16). The main difference is that the hazard in Eq. (5.19) jumps from the state of absolute refractoriness to a constant firing rate, whereas in Eq. (5.16) the transition is smooth.

## 5.2.4 Stationary renewal theory and experiments

Renewal theory is usually associated with stationary input conditions. The interval distribution P0 can then be estimated experimentally from a single long spike train. The applicability of renewal theory relies on the hypothesis that a memory back to the last spike suffices to describe the spike statistics. In particular, there should be no correlation between one interval and the next. In experiments, the renewal hypothesis, can be tested by measuring the correlation between subsequent intervals. Under some experimental conditions, correlations are small indicating that a description of spiking as a stationary renewal process is a good approximation (Goldberg et al., 1964).

The notion of stationary input conditions is a mathematical concept that cannot be easily translated into experiments. With intracellular recordings under in vitro conditions, constant input current can be imposed and thus the renewal hypothesis can be tested directly. Under in vivo conditions, the assumption that the input current to a neuron embedded in a large neural system is constant (or has stationary statistics) is questionable; see (Perkel et al., 1967a,b) for a discussion. While the externally controlled stimulus can be made stationary (e.g., a grating drifting at constant speed), the input to an individual neuron is out of control.

Let us suppose that, for a given experiment, we have checked that the renewal hypothesis holds to a reasonable degree of accuracy. From the experimental interval distribution P0 we can then calculate the survivor function S0 and the hazard via Eqs. (5.5) and  (5.10); see the examples in subsection 5.2.2. If some additional assumptions regarding the nature of the noise are made, the form of the hazard (t|) can be interpreted in terms of neuronal dynamics. In particular, a reduced hazard immediately after a spike is a signature of neuronal refractoriness (Goldberg et al., 1964; Berry and Meister, 1998).

In case of a stationary renewal process, the interval distribution P0 contains all the statistical information, in particular mean firing rate, autocorrelation function and noise spectrum can be derived.

### 5.2.4.1 Mean firing rate

To arrive at an expression for the mean firing rate, we start with the definition of the mean interval,

 s = s P0(s) ds . (5.21)

The mean firing rate has been defined in Chapter 1.4 as = 1/s. Hence,

 = s P0(s) ds  = S0(s) ds . (5.22)

The second equality sign follows from integration by parts using P0(s) = - dS0(s)/ds; cf. Eq. (5.5).

### 5.2.4.2 Autocorrelation function

Let us consider a spike train Si(t) = (t - ti(f)) of length T. The firing times ti(f) might have been measured in an experiment or else generated by a neuron model. We suppose that T is sufficiently long so that we can formally consider the limit T. The autocorrelation function Cii(s) of the spike train is a measure for the probability to find two spikes at a time interval s, i.e.

 Cii(s) = Si(t) Si(t + s) , (5.23)

where . denotes an average over time t,

 f (t) = f (t) dt . (5.24)

We note that the right-hand side of Eq. (5.23) is symmetric so that Cii(- s) = Cii(s) holds.

The calculation of the autocorrelation function for a stationary renewal process is the topic of the next section.

### 5.2.4.3 Noise spectrum

The power spectrum (or power spectral density) of a spike train is defined as () = limTT(), where T is the power of a segment of length T of the spike train,

 T() = Si(t) e-i t dt , (5.25)

The power spectrum () of a spike train is equal to the Fourier transform () of its autocorrelation function (Wiener-Khinchin Theorem). To see this, we use the definition of the autocorrelation function

 () = Si(t) Si(t + s) e-i s ds =  Si(t)  Si(t + s) e-i s ds dt =  Si(t) e+i t dt Si(s') e-i s' ds' . (5.26)

In the limit of T, Eq. (5.25) becomes identical to (5.26) so that the assertion follows. The power spectral density of a spike train during spontaneous activity is called the noise spectrum of the neuron (Bair et al., 1994; Edwards and Wakefield, 1993). As we will see in the next subsection, the noise spectrum of a stationary renewal process is intimately related to the interval distribution P0(s).

## 5.2.5 Autocorrelation of a stationary renewal process

Noise is a limiting factor to all forms of information transmission and in particular to information transmission by neurons. An important concept of the theory of signal transmission is the signal-to-noise ratio. A signal that is transmitted at a certain frequency should be stronger than (or at least of the same order of magnitude as) the noise at the same frequency. For this reason, the noise spectrum () of the transmission channel is of interest. In this section we calculate the noise spectrum of a stationary renewal process. As we have seen above, the noise spectrum of a neuron is directly related to the autocorrelation function of its spike train. Both noise spectrum and autocorrelation function are experimentally accessible (Bair et al., 1994; Edwards and Wakefield, 1993).

Let = Si denote the mean firing rate (expected number of spikes per unit time) of the spike train. Thus the probability of finding a spike in a short segment [t, t + t] of the spike train is  t. For large intervals s, firing at time t + s is independent from whether or not there was a spike at time t. Therefore, the expectation to find a spike at t and another spike at t + s approaches for s a limiting value Si(tSi(t + s) = Cii(s) = . It is convenient to subtract the baseline value and introduce a `normalized' autocorrelation,

 C0ii(s) = Cii(s) -  , (5.27)

with Cii0(s) = 0. Fourier transform of Eq. (5.27) yields

 () = () + 2 () . (5.28)

Thus () diverges at = 0; the divergence is removed by switching to the normalized autocorrelation. In the following we will calculate () for 0.

In the case of a stationary renewal process, the autocorrelation function is closely related to the interval distribution P0(s). This relation will now be derived. Let us suppose that we have found a first spike at t. To calculate the autocorrelation we need the probability density for a spike at t + s. Let us construct an expression for Cii(s) for s > 0. The correlation function for positive s will be denoted by  C+(s) or

 C+(s) =  Cii(s) (s) . (5.29)

The factor in Eq. (5.29) takes care of the fact that we expect a first spike at t with rate . C+(s) gives the conditional probability density that, given a spike at t, we will find another spike at t + s > t. The spike at t + s can be the first spike after t, or the second one, or the nth one; see Fig. 5.4. Thus for s > 0
 C+(s) = P0(s) + P0(s') P0(s - s') ds' + P0(s') P0(s'') P0(s - s' - s'') ds' ds'' + ... (5.30)

or

 C+(s) = P0(s) + P0(s') C+(s - s') ds' (5.31)

as can be seen by inserting Eq. (5.30) on the right-hand side of (5.31).

Due to the symmetry of Cii, we have Cii(s) =  C+(- s) for s < 0. Finally, for s = 0, the autocorrelation has a peak reflecting the trivial autocorrelation of each spike with itself. Hence,

 Cii(s) = (s) + C+(s) + C+(- s) . (5.32)

In order to solve Eq. (5.31) for C+ we take the Fourier transform of Eq. (5.31) and find

 () =  , (5.33)

Together with the Fourier transform of Eq. (5.32), =  [1 + 2 Re{C+()}], we obtain

 () =  Re    for    0 . (5.34)

For = 0, the Fourier integral over the right-hand side of Eq. (5.30) diverges, since P0(s)ds = 1. If we add the diverging term from Eq. (5.28), we arrive at

 () =  Re +2 () (5.35)

This is a standard result of stationary renewal theory (Cox and Lewis, 1966) which has been repeatedly applied to neuronal spike trains (Bair et al., 1994; Edwards and Wakefield, 1993).

### 5.2.5.1 Example: Stationary Poisson process

In Section 5.2.3 we have defined the Poisson neuron as a stationary renewal process with constant hazard (t - ) = . In the literature, a Poisson process is often defined via its autocorrelation

 Cii(s) =  (s) + (5.36)

We want to show that Eq. (5.36) follows from Eq. (5.30).

Since the interval distribution of a Poisson process is exponential [cf. Eq. (5.18)], we can evaluate the integrals on the right-hand side of Eq. (5.30) in a straightforward manner. The result is

 C+(s) =  e- s1 +  s + ( s)2 + ... =  . (5.37)

Hence, with Eq. (5.32), we obtain the autocorrelation function (5.36) of a homogeneous Poisson process. The Fourier transform of Eq. (5.36) yields a flat spectrum with a peak at zero:

 () = +2  () . (5.38)

The result could have also been obtained by evaluating Eq. (5.35).

### 5.2.5.2 Example: Poisson process with absolute refractoriness

We return to the Poisson neuron with absolute refractoriness defined in Eq. (5.19). Apart from an absolute refractory time , the neuron fires with rate r. For 0, Eq. (5.35) yields the autocorrelation function

 () =  1 + 2  [1 - cos( )] + 2 sin( ) , (5.39)

cf. Fig. (5.4)B. In contrast to the stationary Poisson process Eq. (5.36), the noise spectrum of a neuron with absolute refractoriness > 0 is no longer flat. In particular, for 0, the noise level is decreased by a factor [1 + 2( ) + ( )2]-1. Eq. (5.39) and generalizations thereof have been used to fit the power spectrum of, e.g., auditory neurons (Edwards and Wakefield, 1993) and MT neurons (Bair et al., 1994).

Can we understand the decrease in the noise spectrum for 0? The mean interval of a Poisson neuron with absolute refractoriness is s = + r-1. Hence the mean firing rate is

 =  . (5.40)

For = 0 we retrieve the stationary Poisson process Eq. (5.2.3) with = r. For finite the firing is more regular than that of a Poisson process with the same mean rate . We note that for finite > 0, the mean firing rate remains bounded even if r. The neuron fires then regularly with period . Because the spike train of a neuron with refractoriness is more regular than that of a Poisson neuron with the same mean rate, the spike count over a long interval, and hence the spectrum for 0, is less noisy. This means that Poisson neurons with absolute refractoriness can transmit slow signals more reliably than a simple Poisson process.

Next: 5.3 Escape noise Up: 5. Noise in Spiking Previous: 5.1 Spike train variability
Gerstner and Kistler
Spiking Neuron Models. Single Neurons, Populations, Plasticity
Cambridge University Press, 2002