International Journal of Bifurcation and Chaos, Vol. 14, No. 6 (2004) 2061–2068 c World Scientific Publishing Company
RELIABILITY OF SPIKE TIMING IN A NEURON MODEL J. M. CASADO ´ Area de F´ısica Te´ orica. Universidad de Sevilla, Apartado Correos 1065, 41080 Sevilla, Spain ´ J. P. BALTANAS Departamento de Matem´ aticas y F´ısica Aplicadas y Ciencias de la Naturaleza, Escuela Superior de Ciencias Experimentales y Tecnolog´ıa, Universidad Rey Juan Carlos, 28933 M´ ostoles, Spain Received October 29, 2001; Revised April 12, 2002 Recently, it has been shown experimentally [Mainen & Sejnowski, 1995] that, in contrast to the lack of precision in spike timing associated with flat (dc) stimuli, neocortical neurons of rats respond reliably to weak input fluctuations resembling synaptic activity. This has led authors to suggest that, in spite of the high variability of interspike intervals found in cortical activity, the mechanism of spike generation in neocortical neurons has a low level of intrinsic noise. In this work we approach the problem of spike timing by using the well-known FitzHugh–Nagumo (FHN) model of neuronal dynamics and find that here also, fluctuating stimuli allow a more reliable temporal coding than constant suprathreshold signals. This result is associated with the characteristics of a phenomenological stochastic bifurcation taking place in the noisy FHN model. Keywords: Neural computation; neural code; encitable systems; Hopf bifurcation; reliability; spike timing.
1. Introduction Neurons transmit information by encoding timedependent stimuli into sequences of discrete, stereotyped action potentials or spikes. Early theories assumed that all the information about the input signals is carried by the firing rate and that the observed variability in spike timing is just a sort of noise. This is the basic idea behind the concept of rate coding. Nevertheless, in the last few years, more and more studies have shown that irregularity in the timing of individual spikes represents the basic mechanism of long-distance information transmission along the neural pathways [Rieke et al., 1996]. If this is so, the analysis of the temporal precision with which neurons are capable of encoding a stimulus into a spike train is a crucial issue and
deserves careful consideration. Recently, Mainen and Sejnowski [1995] have compared the response of neocortical neurons of rats to different stimuli and have shown that, in contrast with the lack of precision of the firing patterns evoked with flat (dc) current pulses, fluctuating input signals resembling synaptic activity produce highly reliable spike sequences. Their data suggest that, in spite of the high variability present in cortical activity, the noise present within neurons of the neocortex seems to be small, and furthermore, the noise seems to play a much smaller role in the detection of “natural” signals, like fluctuating pulses mimicking synaptic activity, than in the case of flat pulses coming from laboratory sources. Indeed, a very similar result was reported in crayfish neurons two decades earlier by Bryant and Segundo [1976]. 2061
2062
J. M. Casado & J. P. Baltan´ as
Neuronal noise has been usually associated with a variety of different factors, the most important of them probably being the fluctuations in synaptic input due to the activity of adjacent nerve cells. In mathematical models of excitable dynamics, the noise is usually introduced by means of stochastic terms with some predetermined correlation properties. In the paper cited above, Bryant and Segundo [1976] remarked that their results could be explained by a simple linear filter and threshold crossing. Recently, Schneidman et al. [1998] have been able to show that the results of Mainen and Sejnowski [1995] are attributable to fluctuations in the small number of Na+ /K+ channels that are open near threshold. The aim of this paper is to show that the high reproducibility of the spike train induced in a neuron by fluctuating stimuli can be understood by using a much simpler model of spiking dynamics.
2. The Stochastic FitzHugh Nagumo Model The FHN model is described by the following pair of differential equations dv = v(v − v0 )(1 − v) − w + A , dt (1) dw τw = v − w − a0 , dt where a0 and v0 (0 < v0 < 1) are some constants and the parameters τv and τw allow us to separate the characteristic time scales of v(t) and w(t). In the context of neurophysiology, v(t) describes the voltage across the neuron’s membrane, whereas w(t) corresponds to the recovery of the polarity of the membrane after the firing of an action potential. The excitable character of the FHN model relies on the existence of a threshold where a Hopf bifurcation takes place. The precise localization of the bifurcation point depends on the numerical values of the various parameters appearing in Eq. (1). Here we shall restrict ourselves to v0 = 0.5 and a0 = 0.15 and, furthermore, we shall choose the time scales to be τv = 0.01 and τw = 1, so that the voltage variable is much faster than the recovery variable. It is easy to show that for such a set of parameter values, the system has a globally stable fixed point as long as A < Ac ≈ 0.115. For A greater than this value, the fixed point becomes unstable and a stable limit cycle starting from amplitude zero appears. A numerical simulation of Eq. (1) shows τv
Fig. 1. Bifurcation diagrams showing the amplitude and period of oscillations of v(t) as functions of parameter A for the FHN model. For values of A just past the bifurcation point, the system performs low amplitude (subthreshold) oscillations. However, a sufficiently strong constant signal activates a relaxational oscillation of greater amplitude giving rise to the firing of a spike. The continuous lines are plotted only as a guide to the eye.
that for A very close to Ac (with A > Ac ) the initial low-amplitude periodical oscillations abruptly transform into large-amplitude, relaxational excursions around the phase plane. In Fig. 1 the amplitude of the oscillations of the voltage variable past the supercritical Hopf bifurcation is depicted as a function of the parameter A. Recently, Osipov and Ponizovskaya have found an analogous behavior in a two-dimensional approximation of the Hindmarsh– Rose model [Osipov & Ponizovskaya, 1998]. Let us consider now the stochastic version of the FHN model, τv
dv = v(v − v0 )(1 − v) − w + A(t) , dt
dw = v − w − a0 , τw dt
(2)
Reliability of Spike Timing in a Neuron Model
Fig. 2. A short section of the spike train generated by the FHN system under the action of a fluctuating input. Observe the induced switching from low amplitude oscillations to spikes and back. The short vertical lines correspond to the firing times.
where the total input A(t) is supposed to be the sum of two contributions A(t) = S(t) + ξ(t) ,
(3)
S(t) representing a time-dependent stimulus, whereas ξ(t), the stochastic process describing the neural noise, is assumed to be a Gaussian, zero mean, stochastic process with a time correlation function given by hξ(t)ξ(s)i = Dδ(t − s) ,
(4)
where t and s are two arbitrary instants of time and D is the noise intensity. In the experiments performed by Mainen and Sejnowski [1995], an exponentially filtered, Gaussian white noise signal was delivered to the neuron in order to simulate the net current injected by adjacent nerve cells. Sequences of this computer generated signal were delivered after adding it to a constant depolarizing pulse in order to measure the reliability of the evoked response. It is then reasonable, in the context of the somewhat crude model being used here, to describe fluctuations of ionic channels near the threshold of firing as a stochastic process that is added to the injected signal. It has been shown that beyond the supercritical Hopf bifurcation, the response of the FHN system to fluctuating stimuli can display a random pattern of transitions between quasi-harmonic, small-amplitude oscillations and large-amplitude, relaxational excursions of the phase variables
2063
[Makarov et al., 2001]. This dynamical bistability could be associated with the particular form in which the amplitude of the oscillatory motion increases past the bifurcation threshold, as suggested by the work of Osipov and Ponizovskaya [1998]. To the best of our knowledge, however, a detailed analysis of this important problem has not been carried out. Each one of those large relaxational excursions can be associated with the firing of one spike and, thus, it is possible to characterize the behavior of the forced system by analyzing the spike train induced by the fluctuating stimulus. In Fig. 2, a short section of the typical spike train generated by the FHN model is presented to show its characteristic features. In this case, the stimulus was provided by a stochastic process ξ(t) with hξ(t)i < A c .
3. Numerical Results In this section, we consider the stochastic FHN model described above and study the response of the system when it is stimulated by both a flat and a fluctuating signal. In particular, we follow the work of Mainen and Sejnowski [1995] and focus our attention on the characterization of the reproducibility of the firing times of the spikes triggered by the system after the onset of the external stimulus. Let us consider first the case in which we stimulate the system with a flat pulse S(t) that takes the value S1 > Ac for some interval (t0 , t1 ), being S0 < Ac otherwise [see Fig. 3(a)]. Starting from given values of v(0) and w(0), the equations of the model can be solved by using a numerical algorithm, thus obtaining the response of the system to the stimulus. All the results presented in this work have been obtained by using a fourth order, stochastic Runge–Kutta algorithm with a time step of h = 0.001 [Honeycutt, 1992]. As S 0 = −0.2 is well below the threshold for firing, the system variables relax initially to the stable fixed point and remain in its vicinity. When the pulse is switched on by setting S1 = 0.15 at t0 , the system enters the suprathreshold regime and starts an oscillatory motion. In absence of noise (D = 0), the flat pulse triggers a periodic firing activity with a period that is dependent on the value S1 . Upstrokes of the voltage variable are counted as spikes when they reach a threshold at vc = 0.75, having previously crossed a reset value of 0.2 from below. This spike detection scheme discards any rapid recrossing of the
2064
J. M. Casado & J. P. Baltan´ as
0.3
60
-0.2 -0.1 -0.3 -0.2
-0.3 0.3
S(t) S(t)
Firing rate
0 0.1 -0.1 0
0.2 0.3 0.1 0.2 0 0.1 -0.1 0 -0.2 -0.1 -0.3 -0.2
A
40 20 0
10
14
18
22
26
10
14(a) 18
22
26
Firing rate
S(t) S(t)
0.2 0.3 0.1 0.2
B
60 40 20 0 0
10
14
10
14
-0.3
18 Time 18 Time
22
26
22
26
(b) Fig. 3. (a) Flat stimulus, (b) fluctuating stimulus, used in this work to excite the spiking mechanism of the FHN model. The fluctuating signal corresponds to an Ornstein– Uhlenbeck process with µ = 0.15, Ds = 5 × 10−5 and τs = 100h. In both cases, the suprathreshold portion of the pulse starts at t0 = 10 and lasts for 15 units of time. The bifurcation threshold for constant signals is also depicted (dotted line).
threshold and introduces after each one of the firing episodes a refractory period in which the system is unable to be excited. When D 6= 0, there is a noisy component added to the flat pulse and each presentation of the same signal S(t) results in a different realization of the process (v(t), w(t)), thus giving rise to different firing times. In this form, the variability in the spike timing is associated with the existence of noise in our model. The mean number of spikes fired by the system during the time course of the pulse has been obtained from 5000 realizations of the stochastic process v(t), each one following a different presentation of the same stimulus in presence of noise. This allows us to calculate both the mean firing rate of the system and the fluctuations in the number of spikes
2
4
6 8 10 12 14 16 Time
Fig. 4. PSTH associated with the response of the stochastic FHN model to a pulse of current. In both panels D = 7 × 10−7 . (A) A flat pulse of current, (Ds = 0), has been used to activate the spiking mechanism, whereas (B) a fluctuating pulse with Ds = 10−4 has been employed. In both cases 5000 realizations of the pulse were used and the timing of the successive spikes were accumulated in bins of size ∆t = 0.02. To obtain the firing rate, the number of spikes accumulated in each bin was divided by the total number of trials and by the width of the bin. Time is measured starting from t0 . The very sharp peak near zero corresponds to the first spike fired after each presentation of the stimulus.
fired during the suprathreshold section of the stimulus. Also, during each presentation of the stimulus, the successive firing times were determined and its values accumulated in a suitable discretization of the time axis. By counting the number of spikes in each bin and after an adequate normalization, the post-stimulus time histogram (PSTH) can be obtained. In panel A of Fig. 4, such a histogram is depicted to illustrate the response of the FHN model to the aforementioned stimulus. Apart from the first spike of each train, which is tightly locked to the onset of the pulse, the variability of the successive interspike intervals is summed to increase the desynchronization of corresponding action potentials over the time course of the stimulus. This variability comes from variations in the period of the limit cycle past the bifurcation threshold due to the presence of noise. Let us consider now the effect of a fluctuating stimulus. The arrival of many uncorrelated excitatory and inhibitory synaptic inputs to the neuron results in a total current that may be described by
Reliability of Spike Timing in a Neuron Model
0.008
30 55 250 700
25 Firing rate
0.004 Current
2065
0 −0.004
20 15 10 5 0
−0.008 −2
−1.5
−1
−0.5
0
Time
2.6
2.8
3
3.2 Time
3.4
3.6
3.8
Fig. 5. Spike-triggered stimulus averages for different values of τs . Each average was calculated by using 25 different portions of the external stimulus preceding a spike. Notice that the basic characteristic shape of the reverse correlations is maintained for different correlation times, although it is broadened when τs is increased. The reverse correlations or spike-triggered stimulus averages shown here are similar to the first-order response kernel of the neuron (see e.g. [Rieke et al., 1996]).
Fig. 6. Enlarged portion of the PSTH for Ds = 10−3 and D = 10−6 . A threshold set at three times the mean firing rate of the system over the set of trials was used to select the elevations in instantanous firing rate or events (horizontal line). Firings within the temporal limits set by the threshold (dotted vertical lines) were accumulated and the reliability was calculated by dividing this number with that of total spikes fired during the presentation of the whole block of trials. Spikes belonging to the first peak of the PSTH were discarded.
means of a Gaussian process with a finite correlation time. Thus, to model the stimulus, we have used a pulse of the same duration as before but now changing its flat portion by a fluctuating signal of suprathreshold amplitude. In particular we have used the integral algorithm of Fox et al. [1988] to obtain a realization of an Ornstein–Uhlenbeck process of zero mean and correlation function given by Ds −|t|/τs e , (5) C(t) = τs with τs = 100h and different values of Ds . In order to construct the fluctuating stimulus, we have added this realization to the portion of the flat pulse between times t0 and t1 . Thus, during the suprathreshold portion of the pulse, the mean value of the signal S(t) is µ = 0.15. The particular value D s = 0 corresponds to the suprathreshold portion of the signal S(t) being constant. The relevant portion of the stimulus for a representative case is depicted in Fig. 3(b). Following the work of Mainen and Sejnowski [1995], it is worth introducing here some measures of spike timing derived from the PSTH. These authors define an event as the dramatic elevation in the instantaneous firing rate coming from the clustering of spiking episodes at some instants after the presentation of the stimulus. In Fig. 6, two of those events are clearly identified in the PSTH. In order
to quantify the ability of a neuron to reproduce the same spike train after the presentation of the same stimulus, the authors have introduced two measures termed reliability and temporal precision of the response. The last one is defined as the standard deviation of spike times within any event, averaged over all events during a response. As is clearly shown by the two Post Stimulus Time Histograms depicted in Fig. 4, the response to a flat stimulus has a lower precision than the corresponding to a fluctuating one. In this paper, we will concentrate mainly in the calculation of the reliability associated with the response of the FHN model. For fluctuating stimuli also, both the mean number of spikes fired during the presentation of the stimulus and the PSTH were computed for different values of the intensity of the intrinsic noise. In this case, the bifurcation threshold Ac is only a reference, as there is no bifurcation from a fixed point to a limit cycle. The use of a spike-triggered stimulus average or inverse correlation can reveal the shape of the stimulus that tends to precede the firing of a spike and can indicate the length of stimulus that is relevant to the spiking process. To compute spiketriggered averages, we have averaged over a great deal of firing episodes the form of the signal immediately preceding the firing of a spike. As shown in Fig. 5, the firing episodes are associated with rapid increases subsequent to negative departures of S(t)
J. M. Casado & J. P. Baltan´ as
from its mean value, just as obtained by Mainen and Sejnowski [1995] in their experiments. In contrast with the firings triggered by the flat pulse, now almost every spike is the consequence of the noiseinduced switching from the small-amplitude to the large-amplitude regimes. Thus, spiking is governed by the temporal structure of the stimulus, the noise playing a secondary role unless its intensity is comparable to that of the stimulus. Recently, Tanabe and Pakdaman [2001] have shown that the structure of the phenomenological stochastic bifurcation characterizing the transition from oscillatory to relaxational oscillations resembles closely the Hopf bifurcation of the deterministic model shown in Fig. 1. Their results suggest that the stochastic bifurcation of the noisy FHN system is in fact an extension of the Hopf bifurcation that takes place in the system of deterministic equations governing the behavior of the moments in the Gaussian approximation. As the intervals between firing times are quite independent of low-intensity noise, the reproducibility of the spike timing in the case of fluctuating stimulus is higher than that obtained for the flat stimulus where the intervals between firing times were only dependent on noise. As the firings are triggered by changes of regime, if the fluctuating stimulus is too weak to induce those changes, the reliability will be analogous to that obtained for a flat pulse. The length of the portion of the signal relevant to the firing of a spike is associated with the correlation time of the Ornstein–Uhlenbeck process from which the stimulus is drawn. For a realization of white noise, for example, this length will be trivially zero. This result, which is detailed in Fig. 5, also shows qualitative agreement with the findings of Mainen and Sejnowski [1995]. As can be seen in the PSTH depicted in Fig. 4(B), the spikes fired by the system under the action of a fluctuating stimulus show very little variability from trial to trial and the width of the peaks around the deterministic (D = 0) firing times does not increase along the spike train. We have used a threshold set at three times the mean firing rate of the system during all the presentations of one stimulus to select the periods of elevated firing rate or events. The procedure is outlined in Fig. 6. Then, the reliability is defined as the fraction of total spikes taking place during such periods of elevated firing rate. In Fig. 7, we have plotted the reliability as a function of the signal strength and for three
1 0.8 Reliability
2066
0.6
a
0.4
b
0.2
c
0
↑
-6
-5 log Ds
-4
-3
Fig. 7. Reliability of the response of the FHN system to a suprathreshold signal as a function of its strength and for some values of the noise strength. (a) D = 3 × 10−7 , (b) D = 7 × 10−7 , (c) D = 10−6 . To calculate the reliability, the first spike of each train was discarded. The arrow indicates the values corresponding to a flat signal (Ds = 0). Lines joining different symbols are plotted as a guide for the eye.
different noise intensities. In each case, the reproducibility of the spike train increases when the amplitude of the stimulus is increased. On the other hand, we observe that, for a given stimulus, the spike pattern fired by the FHN model increases its reliability as long as the noise is lowered. For D = 3 × 10−7 , there is a slight decrease in reliability as Ds is increased from 0 to 10−5 . This behavior is associated with the particular sample of signal being used. In fact, the particular values of the reliability obtained depend on the samples of the Ornstein–Uhlenbeck process being used as stimulus. In order to get rid of this dependence we have employed five different portions of the signal to average over. The results of this averaging for a noise level of D = 10−6 are shown in Fig. 8. In this plot, points represent the mean reliabilities over the five samples and the lengths of the error bars are equal to twice the standard deviations. Again, the reliability of spike patterns is strongly correlated with the amplitude of stimulus fluctuations. As can be observed in Fig. 9, the mean firing rate, i.e. the number of spikes fired in unit time during the suprathreshold portion of the stimulus, averaged over all the different realizations of the signal, is quite independent of the amplitude of the stimulus, provided that this is of moderate intensity. In this plot, error bars measure the variability of the
Reliability of Spike Timing in a Neuron Model
Anyway, as these and other authors have pointed out, they are consistent with several experimental and theoretical studies on cortical information processing in which spike timing is important [Mainen & Sejnowski, 1995; Rieke et al., 1996].
1
Reliability
0.8 0.6 0.4 0.2 0
4. Conclusions ↑
-6
-5 log D s
-4
-3
Fig. 8. Reliability of the response of the FHN model to suprathreshold signals averaged over different sections of the signal. Here each point has been obtained by averaging the reliability over five different sections of the same signal. The noise strength is D = 10−6 . In each case, the error bars are equal to twice the standard deviation. As before, the arrow corresponds to a flat signal (Ds = 0) and lines joining the open circles are plotted as a guide for the eye.
1.2 r-σ/r
2067
0.8 0.4 0
The comparison of our results with those of Mainen and Sejnowski [1995] suggests that the reliability of spike timing found in neurons of the neocortex seems to be quite independent of the detailed model used to describe the excitable dynamics of the neuron. Also, we have found that the effects of additive noise with regard to the timing of spikes are strongly dependent on the kind of signal that stimulates the neuron. This could be important in studies in which neurons are stimulated by a signal embedded in a noisy environment in order to analyze the effect of noise in the signal transduction efficiency of the neuron. Finally, computation of spike-triggered stimulus averages has revealed that both the typical shape of the signal preceding a spike and its dependence on the correlation time of the external signal are consistent with experimental findings.
Acknowledgments ↑
-6
-5 log Ds
-4
-3
Fig. 9. Average firing rate (open circles) and its relative standard deviation (black circles), as functions of the amplitude of the stimulus. 5000 different presentations of the same signal were employed to determine the statistics of firing episodes during the suprathreshold portion of the stimulus. In addition, an average over five different signals with the same characteristics was carried out. Error bars correspond to twice the standard deviation associated with this last average.
spike count associated with the different samples of signal being used. Comparison of these results with those of the reliability displayed in Fig. 6, indicates that, at least for the regime of amplitudes studied, the timing of spikes produced by the FHN model is far more correlated with the stimulus than the mean firing rates. We would like to emphasize that both the results found by Mainen and Sejnowski [1995] and those described here do not prove the physiological relevance of the timing of individual spikes.
J. M. Casado acknowledges the Direcci´on General de Investigaci´on Cient´ıfica y T´ecnica (DGICYT) of Spain (Project PB98-1120), for support. J. P. Baltan´as acknowledges support from the Ministerio de Ciencia y Tecnolog´ıa of Spain under project BFM2000-0967.
References Bryant, H. L. & Segundo, J. P. [1976] “Spike initiation by transmembrane current: A white noise analysis,” J. Physiol. 260, 279–314. FitzHugh, R. [1961] “Impulses and physiological states in models of nerve membrane,” Biophys. J. 1, 445–466. Fox, R. F., Gatland, I. R., Roy, R. & Vemuri, G. [1988] “Fast, accurate algorithm for numerical simulation of exponentially correlated noise,” Phys. Rev. A38, 5938–5940. Honeycutt, R. L. [1992] “Stochastic Runge–Kutta algorithms II. Colored noise,” Phys. Rev. A45, 604–610. Mainen, Z. F. & Sejnowski, T. J. [1995] “Reliability of spike timing in neocortical neurons,” Science 268, 1503–1506.
2068
J. M. Casado & J. P. Baltan´ as
Makarov, V. A., Nekorkin, V. I. & Velarde, M. G. [2001] “Spiking behavior in a noise-driven system combining oscillatory and excitatory properties,” Phys. Rev. Lett. 86, 3431–3434. Nagumo, J., Arimoto, S. & Yosizawa, Proc. IRE [1962] 50 in FitzHugh, R. [1962] Biological Engineering, ed. Schwan, H. P. (McGraw-Hill, NY), p. 2061. Osipov, V. V. & Ponizovskaya [1998] “The nature of bursting noises, stochastic resonance and deterministic chaos in excitable neurons,” Phys. Lett. A238, 369–374.
Rieke, F., Warland, D., de Ruyter van Steveninck, R. & Bialek, W. [1996] Spikes. Exploring the Neural Code (MIT Press, Cambridge, MA). Schneidman, E., Freedman, B. & Segev, I. [1998] “Ionchannel stochasticity may be critical in determining the reliability and precision of spike timing,” Neural Comput. 10, 1679–1703. Tanabe, S. & Pakdaman, K. [2001] “Dynamics of moments of FitzHugh-Nagumo neuronal models and stochastic bifurcations,” Phys. Rev. E63, 031911.