Analysis of noise-induced temporal correlations in neuronal spike sequences

UPCommons Portal del coneixement obert de la UPC http://upcommons.upc.edu/e-prints Aquesta és una còpia de la versió author’s final draft d'un article publicat a la revista European physical journal: special topics . La publicació final està disponible a Springer a través de http://dx.doi.org/10.1140/epjst/e2016-60024-6 This is a copy of the author 's final draft version of an article published in the journal European physical journal: special topics . The final publication is available at Springer via http://dx.doi.org/10.1140/epjst/e2016-60024-6 Article publicat / Published article: Reinoso, J.A., Torrent, M.C., Masoller, C. (2016) Analysis of noiseinduced temporal correlations in neuronal spike sequences. " European physical journal: special topics ". Vol. 225, issue 13. p.2689-2696. Doi: 10.1140/epjst/e2016-60024-6 EPJ manuscript No. (will be inserted by the editor) Analysis of noise-induced temporal correlations in neuronal spike sequences Jos´e A. Reinoso1,a, M.C. Torrent1,b, and Cristina Masoller1,c Departament de Fisica, Universitat Politecnica de Catalunya, Colom 11, ES-08222 Terrassa, Barcelona, Spain Abstract. We investigate temporal correlations in sequences of noiseinduced neuronal spikes, using a symbolic method of time-series analysis. We focus on the sequence of time-intervals between consecutive spikes (inter-spike-intervals, ISIs). The analysis method, known as ordinal analysis, transforms the ISI sequence into a sequence of ordinal patterns (OPs), which are defined in terms of the relative ordering of consecutive ISIs. The ISI sequences are obtained from extensive simulations of two neuron models (FitzHugh-Nagumo, FHN, and integrate- and-fire, IF), with correlated noise. We find that, as the noise strength increases, temporal order gradually emerges, revealed by the existence of more frequent ordinal patterns in the ISI sequence. While in the FHN model the most frequent OP depends on the noise strength, in the IF model it is independent of the noise strength. In both models, the correlation time of the noise affects the OP probabilities but does not modify the most probable pattern. 1 Introduction Excitable systems are ubiquitous in nature, and in these systems, noise plays a crucial role in enhancing the detection of and response to weak external signals [1,2]. In brain neurons, excitable spikes play a fundamental function in information encoding, processing and transmission [3]. How sensory neurons codify external stimuli in spike sequences is not fully understood, and a lot of efforts are devoted to understand temporal correlations in the spikes of individual neurons. Spike correlations have been investigated by analyzing the sequence of timeintervals between consecutive spikes (inter-spike-intervals, ISIs). A well-known measure for detecting temporal correlations in the ISI sequence, {. . . Ij, Ij+1, Ij+2 ...}, is the serial correlation coefficient, ρn=⟨IjIj+n⟩ − ⟨Ij⟩2 ⟨I2 j⟩− ⟨Ij⟩2.(1) ae-mail: [email protected] be-mail: [email protected] ce-mail: [email protected] 2 Will be inserted by the editor Simple neuron models (such as integrate-and-fire and FitzHugh-Nagumo models) driven by white noise generate uncorrelated ISI sequences. However, correlated noise or time-dependent effects (such as fatigue, adaptation, etc.) generate temporally correlated spikes which reproduce observed ISI statistics and correlations [4–6]. Methods based in symbolic time-series analysis have been shown to be able to provide useful complementary insight. A popular symbolic method is ordinal analysis [7–11], in which a time-series is transformed into a sequence of ordinal patterns (OPs) by taking into account the relative ordering of the values in the time-series. The OP probabilities provide complementary information to that gained from the linear correlation coefficient [12–14]. In Ref. [13] we used ordinal analysis to study ISI correlations induced by the interplay of a weak deterministic input (a subthreshold sinusoidal signal) and Gaussian white noise. We found that, as the amplitude of the periodic signal increased (but it was kept below the firing threshold), ordinal analysis allowed to uncover signatures of determinism in the spike sequence: we found the existence of preferred ordinal patterns, which varied with the strength of the noise and with the period of the input signal. We also found a clear hierarchical structure in the probabilities of the OPs, which exhibited a resonance-like behavior: specific periods and noise levels enhanced temporal ordering in the ISI sequence, maximizing the probabilities of the preferred ordinal patterns. Natural questions that arise from the previous study refer to the role of noise. When the spikes are generated only by an stochastic input (correlated noise), which patterns are more/less probable? How do they depend on the noise strength and correlation time?. The aim of this paper is to investigate how temporal ordering in the ISI sequence emerges, when the spikes are purely induced by noise. With this goal we simulate spike sequences using two well-known neuron models: FitzHugh-Nagumo (FHN) and integrate-and-fire (IF) driven by an Ornstein-Uhlenbeck (OU) noise, and study how the OP probabilities depend on the strength and on the correlation time of the noise. This paper is organized as follows. Section II describes the models used, Sec. III describes the method of ordinal analysis, and Sec. IV presents the results. We consider the influence of the noise strength and also discuss the effect of the correlation time. Section V presents the conclusions. 2 Models 2.1 FitzHugh-Nagumo model The FHN model equations are [15]: ϵdx dt =x−x3 3−y, (2) dy dt =x+a+ζ, (3) where xis the fast variable, yis the slow one, ϵ << 1 and ais a control parameter such that, if |a|>1 the model has a stable node and if |a|<1 a limit circle. ζ represents the noise, which is an Ornstein–Uhlenbeck (OU) process, dζ dt =−λζ +λξ(t),(4) Will be inserted by the editor 3 where ξis a zero-mean Gaussian white noise of intensity D.ζis correlated as < ζtζs>=λDe −(s−t) λ, and thus, λ−1is the correlation time and σ2=λD is the noise intensity. In the simulations we keep a= 1.05 and ϵ= 0.01 constant, and vary σ2and λ−1as control parameters. Figure 1 displays a typical time series, where the spike times, {. . . , tj, tj+1, . . .}, are detected by using a threshold. Then, the sequence of inter-spike-intervals (ISIs) is {. . . Ij=tj+1 −tj. . .}. Fig. 1. (Color online) Time evolution of y(t) obtained from simulation of the FHN model. The circles indicate the detection of the spikes, using as a threshold y= 1.5. The noise intensity is σ2= 0.03 and the correlation time is λ−1= 2. 2.2 Integrate and fire model The sub-threshold dynamics of an IF neuron is governed by [3] ˙v=b−v+ζ, (5) where bis a constant base current and ζis the noise. In this model, when the voltage variable, v, grows above a threshold value, vt, a spike is fired and vis re-set to v0. Without noise, the IF model has a stable fixed point if b < 1 and produces periodic spikes if b > 1. The parameters used in the simulations are: b= 0.97, vt= 1 and v0= 0. ζis a correlated OU-noise, defined as in the FHN model; the noise strength, σ2, and the correlation time, λ−1, are control parameters. Figure 2 displays a typical sequence of spikes. By definition, the spike times, {. . . , tj, tj+1, . . .}, are given by the condition v=vt. 3 Ordinal analysis Ordinal analysis transforms the sequence of time-intervals between consecutive spikes (ISI sequence), {. . . Ij, Ij+1, Ij+2 . . .}into a sequence of symbolic patterns (referred to as ordinal patterns), by taking into account the relative ordering of Lconsecutive values. As schematically represented in the inset of Fig. 3(a), for patterns of length L= 3: 4 Will be inserted by the editor 6220 6230 6240 6250 6260 6270 6280 0 0.5 1 1.5 2 2.5 3 3.5 time v Fig. 2. (Color online) Sequence of spikes simulated with the IF model. The noise intensity is σ2= 0.025 and the correlation time is λ−1= 2. (Ij< Ij+1 < Ij+2) gives 012, label i= 1, (Ij+1 < Ij< Ij+2) gives 102, label i= 2, . . . (Ij+2 < Ij+1 < Ij) gives 210, label i= 6. If two intervals are equal, a tiny random value is added before computing the pattern. This symbolic transformation keeps the information about the temporal ordering in the ISI sequence, but neglects the information contained in the duration of the intervals between consecutive spikes. Next, the OP probabilities are estimated by computing the frequency of occurrence of the different patterns. In order to compute these probabilities with accurate precision, very long time series were simulated, to obtain ISI time series with N= 100000 data points; the influence of the length, N, is discussed at the end of Sec. 4.2. To capture the existence of underlying temporal ordering in the ISI sequence, a significance analysis of the OP probabilities, piwith i= 1 . . . 6, is needed: they have to be not consistent with the uniform distribution. Thus, the null hypotheses (N.H.) is equally probable OPs: pi= 1/6 for i= 1 ...6. To estimate the interval of probability values which is consistent with the N.H., we use a binomial test [16]: considering a confidence level of 95%, if all the OP probabilities are within the range, p±3σp, where p= 1/6 and σp=√p(1 −p)/N, the OP probabilities are consistent with the N.H. (equally probable OPs); thus, there is no preferred temporal ordering in the ISI sequence. In contrast, if at least one probability value is above p+ 3σpor below p−3σp, we infer (with 95% confidence level) that the OP probabilities are not uniformly distributed, i.e., there is significant evidence of underlying preferred (or unfavored) temporal ordering in the ISI sequence. 4 Results The FHN and IF models were simulated using a second-order Runge-Kutta method, with integration step 0.005 for the FHN model and 0.01 for the IF model. The simulations were long enough to generate time-series with 100.000 ISIs. 4.1 Analysis of ordinal ISI correlations in the FitzHugh-Nagumo model Figure 3 displays the OP probabilities vs the noise strength, computed from ISI sequences simulated with the FHN model. Three noise correlation times are considered. Will be inserted by the editor 5 Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.16 0.17 0.18 0.19 0.2 Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.16 0.165 0.17 0.175 Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.16 0.165 0.17 0.175 Fig. 3. (Color online) Probabilities of the ordinal patterns vs the noise strength, σ2, for three values of the correlation parameter, λ: in (a), λ= 0.5; in (b), λ= 1; and in (c), λ= 1.5. In panel (a) the inset indicates the six possible OPs of length L= 3. In each panel, the gray area indicates the region in which the probabilities are consistent with the uniform distribution (pi= 1/6∀i) with 95% confidence level. It can be observed that, for long correlation time (panel a) the OP probabilities are outside the gray region and reveal the existence of two favored OPs: 012 (blue) and 210 (red). The most probable OP depends on the noise strength: is 012 for weak noise and is 210 for stronger noise. As the correlation time decreases (panel b) the OP probabilities are almost all in the grey region and there are no clearly preferred OPs. If the correlation time is further decreased (panel c) all OP probabilities are in the gray region, consistent with equally probable OPs. In each panel, the gray region indicates the range of probability values that are consistent with the N.H.: if all the OP probabilities are within this range, then, there is no evidence of temporal ordering in the ISI sequence. We see three different scenarios: for long correlation time (panel a), the OP probabilities are outside the gray region and the more frequent OPs are 012 and 210. This indicates that there is significant evidence of preferred temporal ordering in the ISI sequence. Remarkably, the most frequent OP depends on the strength of the noise: it is 012 for weak noise and 210 for stronger noise. As the correlation time decreases (panel b), almost all OP probabilities are in the grey region, indicating that there are no clearly preferred OPs. If the correlation time is further decreased (panel c) all OP probabilities are in the gray region, consistent with a uniform distribution of OP probabilities (i.e., there is no evidence of preferred temporal ordering in the ISI sequence). In this case the noise correlations are too weak to induce temporal ordering and the OPs probabilities resemble those with white noise [13]. Figure 4 displays the OP probabilities vs the noise correlation time, now keeping constant the noise strength. We select three values of the noise strength such that the most frequent OP is different: as can be seen in Fig. 3(a), for σ2= 0.01, 012 is the most frequent OP; for σ2= 0.03, 210 is the most frequent OP; and for σ2= 0.02 both, 012 and 210, are about equally probable. We see in Fig. 4 that the noise correlation time does not change the preferred OP, but it changes the OP probabilities: as λ−1increases, the OPs probabilities gradually emerge from the grey region, and the preferred OPs increase their probabilities while the infrequent OPs decrease them: varying λ−1does not induce any change in the preferred OPs, but only in their probabilities. 4.2 Analysis of ordinal ISI correlations in the IF model Next we present the results of the ordinal analysis when the ISI sequences are simulated with the IF model. Figures 5 and 6 display the OP probabilities when either the noise strength (Fig. 5) or the correlation time (Fig. 6) are varied. In Fig. 5(a) we note that the most frequent OP is 012, and it does not change with the strength of the noise. For shorter correlation time, Fig. 5 (b), 012 remains the 6 Will be inserted by the editor λ−1 OP probabilities 1 2 3 4 5 0.15 0.16 0.17 0.18 0.19 0.2 λ−1 OP probabilities 1 2 3 4 5 0.15 0.16 0.17 0.18 0.19 0.2 λ−1 OP probabilities 1 2 3 4 5 0.15 0.16 0.17 0.18 0.19 0.2 Fig. 4. (Color online) As Fig. 3, but the OP probabilities are plotted vs the noise correlation time, λ−1, keeping constant the noise intensity: in (a) σ2= 0.01, in (b) σ2= 0.02 and in (c) σ2= 0.03. We note that the noise correlation time does not modify which are the more frequent and less frequent ordinal patterns, but varies their probabilities. most probable OP (except for very weak noise). For even shorter correlation time, Fig. 5(c), the six probabilities are within the gray region, consistent with equally probable OPs. In the three panels of Fig. 6 we note that the noise correlation time affects the OP probabilities but it does not change which are the most frequent and the less frequent patterns. Taken together, these results suggest that in the IF model, only positive ordinal correlations arise in the ISI sequence (i.e., the favored OP is almost always 012). We note that the OP probabilities have large variability with σ2and λ, as compared with the FHN model. We interpret this as due to the fact that in the IF model the dynamics is less deterministic. The IF model has only one variable which is governed by a linear stochastic rate-equation that describes the subthreshold dynamics, and lacks an equation governing the supra-threshold dynamics (i.e., the evolution during the spikes). In the IF model, action potentials are not generated explicitly, but the membrane voltage is abruptly decreased after each threshold crossing. In contrast, the FHN model consists of two coupled nonlinear stochastic rate-equations that determine the evolution of two variables (a fast and a slow one) which govern the subthreshold dynamics, and also, the dynamics during and after the spikes. Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.16 0.165 0.17 0.175 0.18 Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.16 0.165 0.17 0.175 Noise strength (σ2) OP probabilities 0.01 0.02 0.03 0.04 0.05 0.162 0.164 0.166 0.168 0.17 0.172 Fig. 5. (Color online) As Fig. 3, but the ISI sequence is simulated with the IF model. The noise correlation parameter is (a) λ= 0.05, (b) λ= 0.1, and (c) λ= 0.5. In panel (a) we note that the most frequent OP is 012, and does not change with the strength of the noise. For shorter correlation time (panel b), 012 remains the most probable OP (except for very weak noise). For even shorter correlation time (panel c) all OP probabilities are within the gray region, consistent with equally probable OPs. We conclude by discussing the influence of the length of the time series. We have previously shown, in simulations of the FHN model, that very long time series are needed for a precise calculation of the OP probabilities [13]. This also occurs in Will be inserted by the editor 7 λ−1 OP probabilities 10 20 30 40 50 0.16 0.165 0.17 0.175 0.18 λ−1 OP probabilities 10 20 30 40 50 0.16 0.165 0.17 0.175 0.18 λ−1 OP probabilities 10 20 30 40 50 0.16 0.165 0.17 0.175 0.18 Fig. 6. (Color online) OP probabilities vs noise correlation time, λ−1, computed from ISI sequences generated with the IF model. The noise intensity is: (a) σ2= 0.01, (b) σ2= 0.02 and (c) σ2= 0.03. As in the FHN model, we note here that the noise correlation time affects the OP probabilities but does not change which are the most frequent and less frequent patterns. simulations of the IF model, as shown in Fig. 7 that displays the OP probabilities vs. the length of the time-series. Length of the time series OP probabilities 104105106 0.155 0.16 0.165 0.17 0.175 0.18 Fig. 7. (Color online) OP probabilities, computed from ISI sequences generated with the IF model, vs the length of the time series (note the logarithmic horizontal axis). The noise parameters are σ2= 0.01 and λ= 0.05 . 5 Conclusions To sum up, we have used ordinal analysis to investigate temporal ordering in sequences of spikes induced by correlated noise. We simulated two popular neuron models: FitzHugh-Nagumo and integrate-and-fire and analyzed their output sequences of inter-spike-intervals (ISIs). We found that both, the noise strength and the correlation time affect the probabilities of the ordinal patterns, as preferred patterns (012 and/or 210) emerge in the ISI sequence, indicating underlying temporal ordering in the timing of the spikes. Patterns 012 and 210 are more probable in simulations of the FHN model. This preference is independent of the correlation time and occurs for various noise strengths (that vary the mean ISI), as shown in Fig. 3(a) and Figs. 4(a), 4(b) and 4(c). Pattern 012 is preferred over 210 in ISI sequences generated with the IF model; this can be clearly seen in Figs. 5 and 6. In simulations of the FHN model there is also a difference in the probabilities of the patterns, but the difference is small [in 8 Will be inserted by the editor Fig. 4(a) pattern 012 is slightly more probable than 210, while the opposite occurs in Fig. 4(c)]. In both models these probability differences are significant because the probabilities have been computed from long time-series of ISIs. The difference of the probabilities of 012 and 210 captures time irreversibility in the dynamics, because, when the ISI time-series is analyzed starting from the end towards the beginning (i.e., {. . . Ij+2, Ij+1, Ij. . .}) pattern 012 transforms into 210 and vice-versa. Thus, if there is time-symmetry, P(012) = P(210), while if this symmetry is broken, the difference |P(012) −P(210)|can be a measure of time irreversibility. Due to the fact that sensory neurons encode external stimuli in sequences of correlated spikes, our results could be relevant for understanding the response of individual neurons to time-varying signals, and how they are encoded and transmitted in correlated sequences of spikes. 6 Acknowledgments This work has been supported in part by the Spanish MINECO (FIS2015-66503- C3-2-P) and ITN NETT (FP7-PEOPLE-2011-ITN-289146). C. M. also acknowledges partial support from ICREA ACADEMIA. References 1. L. Gammaitoni, P. H¨anggi, P. Jung, and F. Marchesoni Rev. Mod. Phys. 70, 223 (1998). 2. B. Lindner, J. Garcia-Ojalvo, A. Neiman, and L. Schimansky-Geier, Phys. Rep. 392, 321 (2004). 3. E. Izhikevich, Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. (The MIT press, Cambridge, MA, 2007). 4. M. J. Chacron, A. Longtin, M. St-Hilaire, and L. Maler, Phys. Rev. Lett. 85, 1576 (2000). 5. J. W. Middleton, M. J. Chacron, B. Lindner, and A. Longtin, Phys. Rev. E 68, 021920 (2003). 6. B. Lindner, Phys. Rev. E. 69, 022901 (2004). 7. C. Bandt and B. Pompe, Phys. Rev. Lett. 88, 174102 (2002). 8. Y. Cao, W. W. Tung, J. B. Gao, V. A. Protopopescu, and L. M. Hively, Phys. Rev. E 70, 046217 (2004). 9. U. Parlitz, S. Berg, S. Luther, A. Schirdewan, J. Kurths, and N. Wessel, Comp. Biol. Med. 42, 319 (2012). 10. M. Zanin, L. Zunino, O. A. Rosso and D.Papo, Entropy 14, 1554 (2012). 11. J. M. Amig´o, K. Keller, and J. Kurths, Eur. Phys. J. Spec. Top. 222, 2 (2013). 12. A. Aragoneses, L. Carpi, N. Tarasov, D. V. Churkin, M. C. Torrent, C. Masoller, and S. K. Turitsyn, Phys. Rev. Lett. 116, 033902 (2016). 13. J. A. Reinoso, M. C. Torrent, C. Masoller arXiv:1510.09035 (2015). 14. C. Bandt arXiv:1411.3904 (2014). 15. A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775 (1997). 16. A. Aragoneses, N. Rubido, J. Tiana-Alsina, M. C. Torrent and C. Masoller, Sci. Rep. 3, 1778 (2013).