Neuronal coupling benefits the encoding of weak periodic signals in symbolic spike patterns

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 Communications in nonlinear science and numerical simulation. http://hdl.handle.net/2117/367392 Article publicat / Published paper: Masoliver, M.; Masoller, C. Neuronal coupling benefits the encoding of weak periodic signals in symbolic spike patterns. "Communications in nonlinear science and numerical simulation", 1 Març 2020, vol. 82, p. 105023:1-105023:9. DOI: <https://doi.org/10.1016/j.cnsns.2019.105023>. © <2020>. Aquesta versió està disponible sota la llicència CC-BY- NCND 4.0 http://creativecommons.org/licenses/by-nc-nd/4.0/ Neuronal coupling benefits the encoding of weak periodic signals in symbolic spike patterns Maria Masoliver1, Cristina Masoller1,∗ Abstract The biophysical mechanisms by which an input signal elicits a neuronal response are well known (sufficiently large inputs change the membrane potential of the neuron and generate electrical pulses, known as action potentials or spikes), yet, a good understanding of how neurons use these spikes to encode the signal information remains elusive. Recent theoretical studies have focused on how neurons encode a weak periodic signal (that by itself is unable to generate spikes) in a noisy environment, where stochastic electrical fluctuations that do not encode any information occur. Analyzing spike sequences generated by individual neurons and by two coupled neurons (that were simulated with the stochastic FitzHugh-Nagumo model), it has been found that the relative timing of the spikes can encode the signal information. Using a symbolic method to analyze the spike sequence, preferred and infrequent spike patterns were detected, whose probabilities vary with both, the amplitude and the frequency of the signal. To investigate if this encoding mechanism is plausible also for neuronal ensembles, here we analyze the activity of a group of neurons, when they all perceive a weak periodic signal. We find that, as in the case of one or two coupled neurons, the probabilities of the spike patterns, now computed from the spike sequences of all the neurons, depend on the signal’s amplitude and period, and thus, the patterns’ probabilities encode the information of the ∗Corresponding author Email addresses: [email protected] (Maria Masoliver), [email protected] (Cristina Masoller) 1Universitat Polit`ecnica de Catalunya. Departament de F´ısica. Edifici Gaia, Rambla de Sant Nebridi 22, 08222, Terrassa, Barcelona, Spain. Preprint submitted to Journal of L A T E X Templates May 7, 2019 arXiv:1905.01933v1 [q-bio.NC] 6 May 2019 signal. We also find that the resonances with the period of the signal or with the noise level are more pronounced when a group of neurons perceive the signal, in comparison with when only one or two coupled neurons perceive it. Neuronal coupling is beneficial for signal encoding as a group of neurons is able to encode a small-amplitude signal, which could not be encoded when it is perceived by just one or two coupled neurons. Interestingly, we find that for a group of neurons, just a few connections with one another can significantly improve the encoding of small-amplitude signals. Our findings indicate that information encoding in preferred and infrequent spike patterns is a plausible mechanism that can be employed by neuronal populations to encode weak periodic inputs, exploiting the presence of neural noise. Keywords: neural coding, excitability, spike train variability, neuronal noise, FitzHugh-Nagumo model, time series analysis, symbolic analysis, ordinal analysis 1. Introduction A mechanical input such as tapping someone’s knee elicits a stretch reflex as a response. The biophysical mechanism is known, the muscle stretches as a consequence of the tapping to the tendon, which triggers the generation of spikes by a sensory neuron, which in turn triggers the generation of spikes by a motor neuron, leading to muscle contraction and causing the lower leg to bounce back [1, 2, 3]. On the other hand, the signal encoding mechanism is also known, the frequency at which the sensory neuron fires encodes the information about how fast the muscle is stretching [4]. In turn, the firing rate of the motor neuron encodes the information about the muscle force when it contracts [5]. This is an example of a neural circuit (two neurons interconnected by a synapse), which uses the firing rate code as a coding scheme for an external input. Yet, neurons encode information of different types of signals using different encoding mechanisms, which are not yet fully understood. Neurons can represent external or internal inputs in the timing of the individual spikes, in the relative timing 2 of the spikes of two or more neurons, in the individual firing rate, in the average firing rate of a population of neurons, in the spike arrival times, among other coding schemes [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 3, 17, 18]. Understanding the neural code is crucial, not only to gain knowledge of the operation of the central nervous system, but also, to advance artificial intelligence systems based on neural networks that use spike-processing operations for classification, pattern recognition, logic operations, etc. [19, 20, 21, 22]. Efforts have focused on understanding the role, on neural coding, of neural noise (stochastic electrical fluctuations that do not encode any information [23]) and spike temporal correlations, in particular, for encoding and processing weak sensory signals. By analyzing the coefficient of variation of the inter-spike interval (ISI) distribution (the standard deviation of the ISI distribution divided by the mean), the well-known phenomena of stochastic resonance and coherence resonance have been found. While stochastic resonance [24, 25] refers to the enhancement of weak signal detection, coherence resonance [26, 27, 28, 29, 30, 31] refers to the regularization of the spike train, for an optimal level of noise. On the other hand, ISI correlations lasting several ISIs have been studied by using the lagged serial correlation coefficient, which measures linear relations between sequential interspike intervals [32, 33, 34, 35, 36, 37]. An alternative technique, known as ordinal analysis [38, 39] has also been used to detect nonlinear ISI correlations. In general terms, ordinal analysis transforms a time series into a sequence of symbols, known as ordinal patterns, considering the temporal order relations among the data points in the time series. A main advantage of the ordinal symbolic approach is that it provides a straightforward way to quantify how much information is contained in a sequence of spikes (i.e., to apply Information Theory to the study of the neural code [40, 41]): by counting the number of times each ordinal pattern appears in the spike sequence, the probabilities of the different patterns can be estimated, providing a quantification of the information content. Ordinal analysis has been widely used to investigate biomedical signals, for example, to quantify directionality of coupling in cardio-respiratory data [42], to characterize neuronal spike 3 trains [43, 44, 45], to distinguish healthy subjects from patients suffering from congestive heart failure [46], to classify neurophysiological data [47, 48, 49, 50], etc. In order to understand how a single neuron can encode a weak periodic input signal (that by itself is unable to generate spikes) in the presence of neural noise, Aparicio Reinoso et al. [51] have applied ordinal analysis to spike sequences simulated with the stochastic FitzHugh-Nagumo (FHN) model [52, 53]. It was found that the periodic signal induces temporal order in the form of more and less expressed ordinal patterns. The probabilities of the patterns encode the signal information, as they depend on both, the amplitude and the period of the signal. In a follow up study [54], the role of a second neuron that does not perceive the weak signal was analyzed. It was found that the signal is still encoded in the form of preferred and infrequent ordinal patterns, whose probabilities again depend on the period and amplitude of the signal. An open question is whether this encoding mechanism can be employed by a population of neurons. To answer this question, here we use the stochastic FHN model to simulate the activity of a group of neurons, when they all perceive a periodic signal that is weak enough such that by itself (in the absence of noise) it is unable to generate spikes. Thus, as in previous studies, the neuronal ensemble encodes the signal in spikes sequences which are generated due to the interplay of the signal and the noise. Our main findings can be summarized as follows: (i) the ensemble is able to encode lower amplitude signals, in comparison with the signal amplitude that can be encoded by a single neuron or by two coupled neurons; (ii) the noise-induced and period-induced resonances (some ordinal patterns probabilities are minimum or maximum for particular values of the noise strength or signal period) observed in one [51] or two coupled neuron [54] become more pronounced for the neuronal ensemble and (iii) just a few connections among the neurons can significantly improve the signal encoding. This paper is organized as follows. Section 2 presents the model equations, Sec. 3 presents the ordinal analysis method, Sec. 4 presents the results and Sec. 5 presents the conclusions. 4 2. Model The FitzHugh-Nagumo model is one of the simplest (and yet quite realistic) models that describe excitable systems [52, 53, 55]. The equations describing the dynamics of an ensemble of coupled neurons are: ˙ui=ui−u3 i 3−vi+a0cos(2πt/T ) + σ ki N X j aij(uj−ui) + √2Dξi(t), i 6=j ˙vi=ui+a. (1) Here Nrefers to the number of neurons, vis known as the inhibitor variable and uis known as the activator variable that represents the evolution of the membrane potential: in the excitable regime, if there is no external perturbation or it is not strong enough to overcome the threshold, the membrane potential is held at the resting potential (i.e., stable fixed point) whereas when there is an external perturbation strong enough to overcome the threshold, the membrane potential performs a spike (i.e., an action potential). Typical parameters in the excitable regime are a= 1.05 and = 0.01. The parameters a0and Trepresent the amplitude and period of an external sinusoidal input, and are chosen such that the signal is sub-threshold: without noise the neurons do not fire spikes. Dξi(t) represents an stochastic term of strength D, which is taken as Gaussian distributed, uncorrelated temporally and across the neuronal ensemble: hξi(t)ξj(t0)i=δij δ(t−t0) with hξi(t)i= 0 and hξ2 i(t)i= 1. The neurons are mutually coupled with gap-junction connections, characterized by symmetric links (aij =aji = 1 if neurons iand jare connected, else aij =aji = 0). The coupling strength of each link is σ; to keep the total coupling strength uniform for all neurons, it is normalized by number of connections, ki=Pjaij. Regarding the coupling topology, we focus on all-to-all coupling (in this case ki=N−1 for all i), but we also consider random connections. This allows us to analyze the influence of the number of links, as neurons iand jare connected with probability pthat is varied between 0 and 1. It will 5 be interesting, for future work, to investigate more realistic topologies with, for example, modular or hierarchical structures. The model equations are simulated, from random initial conditions, using the Euler-Maruyama method with an integration step of dt = 10−3. For each set of parameters, the voltage-like variable of each neuron uiis analyzed and the sequence of inter-spike-intervals (ISIs) is computed, {Iji;Iji= (tj+1 −tj)i} with tjdefined by the condition ui(tj) = 0, considering only the ascensions. 3. Ordinal analysis The ordinal method [38] is used to analyze each ISI sequence. From the sequence {I1,...Ii,...IN}(for clarity, the subindex that labels the neuron is removed) symbols known as ordinal patterns are obtained by comparing D consecutive ISIs, based on their temporal relation. For example, if we set D= 2, the total number of possible ordinal patterns is two: 01, for I1< I2and 10, for I1> I2, while if we set D= 3, we have 3! = 6 possible ordinal patterns: 012 (I3> I2> I1), 021 (I2> I3> I1), 102 (I3> I1> I2), 120 (I2> I1> I3), 201 (I1> I3> I2) and 210 (I1> I2> I3). The number of possible ordinal patterns (i.e., the number of possible temporal relations) is determined by the number of permutations, D!. Using the function defined in [46] the sequence of ordinal patterns is computed. In order to determine if there are some preferred/infrequent patterns in the ISI sequences, ordinal patterns probabilities are calculated, taking together all the ISI sequences. The ordinal probabilities are estimated as pi= Ci/M, where Cirefers to the number of times the i−th pattern appears and M=PD! i=1 Ciis the total number of ordinal patterns. If ordinal patterns are equi-probable it does not exist a preferred order relation among the timing of spikes. Yet, if there are preferred/infrequent ordinal patterns, a non-uniform probability distribution is obtained. In order to distinguish between these two cases (uniform vs. non-uniform ordinal distribution) a binomial test is used: if all the ordinal patterns are within the interval [p−3σp, p + 3σp] with p= 1/D! 6 and σp=pp(1 −p)/M the ordinal probabilities are consistent with the uniform distribution with 99.74% confidence level; else, some patterns are over or less expressed than others, and there is some degree of temporal order in the timing of the spikes. A large number of spikes are needed to precisely estimate the ordinal probabilities (the data requirements for a single neuron were analyzed in [51], see Fig. 6). As long simulations are computationally demanding, here we limit to consider ensembles of up to 50 neurons. We have analyzed the role of the number of neurons, and we expect that our findings will hold for larger ensembles. The simulations are done for a time long enough to obtain a total number of 105spikes. As in [51, 54], we use D= 3. This choice is motivated by the fact that only short ISI correlations are expected since the spikes are noise-induced (the signal by itself does not induce spikes). 4. Results The neuronal ensemble displays different dynamical regimes, depending on the coupling strengh, the signal amplitud and period, the noise strength, and the coupling topology. Figures 1 and 2 display several examples of the dynamics of a group of 50 neurons under different conditions: Fig. 1 shows the activity of an individual neuron (the voltage-like variable of neuron 1), while Fig. 2 displays the raster plot of the ensemble. In panels 1(a) and 2(a) the neurons are uncoupled and no signal is applied, therefore, random spiking activity occurs due to the noise. In panels 1(b) and 2(b) the neurons are mutually coupled, still no signal is applied. Now we see synchronized spiking activity superposed with random spikes. When the periodic signal is applied, we see in panels 1(c)- (f) and 2(c)-(f) that the neurons either fire regular and synchronized spikes, or there is more irregular firing, depending on the period of the signal. In the following we analyze the influence of the different parameters. To stress the role of the number of neurons, we compare the results obtained for 50 neurons with those obtained for only two coupled neurons. 7 Fig. 1: Spiking activity of a neuron when no signal is applied (a0= 0) and the neuron (a) is uncoupled (σ= 0), (b) is coupled to a group of 50 neurons (all to all coupling, σ= 0.05). Activity of the neuron when it is coupled and a sinusoidal signal of amplitude a0= 0.1 and period (c) T= 10, (d) T= 20, (e) T= 40 is applied. The noise level is D= 2.5·10−6. Fig. 2: Raster plots displaying the spiking activity of the group of 50 neurons for the same parameters as in Fig. 1 8 is perceived by just one neuron (or by a subset of neurons) propagates on the whole ensemble, which would give information of how the signal is transmitted. As well, we intent to study how an ensemble of neurons may encode two weak signals. A recent experimental study [17] of how neurons encode simultaneous auditory stimuli has found that some neurons fluctuate between firing rates observed for each individual sound. It would be interesting to compare with our synthetic model, to contribute to advance the understanding of how neuronal systems process information of multiple simultaneous stimuli. 6. Acknowledgments This work was supported in part by Spanish MINECO/FEDER grant FIS2015- 66503-C3-2-P141 and ICREA ACADEMIA, Generalitat de Catalunya. 15 Fig. 6: Probabilities of the ordinal patterns as a function of the number of neurons, N(a, b), of the percentage of links (c, d), and of the coupling strength (e, f) for a0= 0.05 and a0= 0.1, respectively. In panels (a, b, e, f) the neurons are all-to-all coupled, in panels (c, d) the coupling topology is random (starting from uncoupled neurons, links are randomly added until the neurons are all-to-all coupled). In panels (c, d, e, f) N= 50, in all the panels: D= 2.5·10−6and T= 10. 16 Fig. 7: Raster plots displaying the spiking activity of the group of 50 neurons all-to-all coupled over time for (a) σ= 0, (b) σ= 0.01, (c) σ= 0.015 and (d) σ= 0.03. In all the panels: D= 2.5·10−6,T= 10 and a0= 0.1. 17 References [1] K. Pearson, J. Gordon, Spinal reflexes, Principles of Neural Science (2000) 713 – 736. [2] J. Nolte, The human brain: Introduction to its functional anatomy, Mosby, St. Louis (2002). [3] J. J. Knierim, Chapter 19 - information processing in neural networks, in: J. H. Byrne, R. Heidelberger, M. N. Waxham (Eds.), From Molecules to Networks, third edition ed., Academic Press, Boston, 2014, pp. 563 – 589. doi:https://doi.org/10.1016/B978-0-12-397179-1.00019-1. [4] C. C. Hunt, Mammalian muscle spindle: peripheral mechanisms, Physiological Rev. 70 (1990) 643–663. [5] A. W. Monster, H. Chan, Isometric force production by motor units of extensor digitorum communis muscle in man, J. of Neurophysiol. 40 (1977) 1432–1443. [6] B. W. Knight, Dynamics of encoding in a population of neurons, J. Gen. Physiol. 59 (1972) 734766. [7] M. Carandini, F. Mechler, C. S. Leonard, J. A. Movshon, Spike train encoding by regular-spiking cells of the visual cortex, J. Neurophysiol. 76 (1996) 3425–3441. [8] Y. Sakurai, How do cell assemblies encode information in the brain?, Neuroscience & Biobehavioral Reviews 23 (1999) 785 – 796. [9] N. Masuda, K. Aihara, Bridging rate coding and temporal spike coding by effect of noise, Phys. Rev. Lett. 88 (2002) 248101. [10] E. Arabzadeh, S. Panzeri, M. E. Diamond, Whisker vibration information carried by rat barrel cortex neurons, J. Neurosci. 24 (2004) 6011–6020. 18 [11] I. Nelken, G. Chechik, T. D. Mrsic-Flogel, A. J. King, J. W. H. Schnupp, Encoding stimulus information by spike numbers and mean response time in primary auditory cortex, J. Comput. Neurosci. 19 (2005) 199–221. [12] B. B. Averbeck, P. E. Latham, A. Pouget, Neural correlations, population coding and computation, Nat. Rev. Neurosci. 7 (2006) 358366. [13] R. Quian Quiroga, S. Panzeri, Extracting information from neuronal populations: information theory and decoding approaches, Nat. Rev. Neurosci. 10 (2001) 173. [14] X.-J. Wang, Neurophysiological and computational principles of cortical rhythms in cognition, Physiol. Rev. 90 (2010) 1195–1268. [15] A. Kumar, S. Rotter, A. Aertsen, Spiking activity propagation in neuronal networks: reconciling different perspectives on neural coding, Nat. Rev. Neurosci. 11 (2010) 615–627. [16] T. Tchumatchenko, A. Malyshev, F. Wolf, M. Volgushev, Ultrafast population encoding by cortical neurons, J. Neurosci. 31 (2011) 12171–12179. [17] V. C. Caruso, J. T. Mohl, C. Glynn, J. Lee, S. M. Willett, A. Zaman, A. F. Ebihara, R. Estrada, W. A. Freiwald, S. T. Tokdar, J. M. Groh, Single neurons may encode simultaneous stimuli by switching between activity patterns, Nat. Comm. 9 (2018) 2715. [18] E. Lazarov, M. Dannemeyer, B. Feulner, J. Enderlein, M. J. Gutnick, F. Wolf, A. Neef, An axon initial segment is required for temporal precision in action potential encoding by neuronal populations, Sci. Adv. 4 (2018). [19] P. R. Prucnal, B. J. Shastri, T. F. de Lima, M. A. Nahmias, A. N. Tait, Recent progress in semiconductor excitable lasers for photonic spike processing, Adv. Opt. Photon. 8 (2016) 228–299. 19 [20] B. J. Shastri, M. A. Nahmias, A. N. Tait, A. W. Rodriguez, B. Wu, P. R. Prucnal, Spike processing with a graphene excitable laser, Sci. Rep. 6 (2016) 19126. [21] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund, M. Soljai, Deep learning with coherent nanophotonic circuits, Nat. Phot. 11 (2017) 441. [22] P. R. Prucnal and B. J. Shastri, Neuromorphic Photonics, CRC Press, 2017. [23] M. D. McDonnell, L. M. Ward, The benefits of noise in neural systems: Bridging theory and experiment, Nat. Rev. Neurosci. 12 (2011) 415–426. [24] J. F. Lindner, B. K. Meadows, W. L. Ditto, M. E. Inchiosa, A. R. Bulsara, Array enhanced stochastic resonance and spatiotemporal synchronization, Phys. Rev. Lett. 75 (1995) 3–6. [25] E. Yilmaz, M. Uzuntarla, M. Ozer, M. Perc, Stochastic resonance in hybrid scale-free neuronal networks, Physica A 392 (2013) 5735 – 5741. [26] A. S. Pikovsky, J. Kurths, Coherence resonance in a noise-driven excitable system, Phys. Rev. Lett. 78 (1997) 775–778. [27] C. Zhou, J. Kurths, B. Hu, Array-enhanced coherence resonance: Nontrivial effects of heterogeneity and spatial independence of noise, Phys. Rev. Lett. 87 (2001) 098101. [28] O. Kwon, H.-T. Moon, Coherence resonance in small-world networks of excitable cells, Phys. Lett. A 298 (2002) 319 – 324. [29] O. Kwon, H.-H. Jo, H.-T. Moon, Effect of spatially correlated noise on coherence resonance in a network of excitable cells, Phys. Rev. E 72 (2005) 066121. [30] P. Balenzuela, P. Ru´e, S. Boccaletti, J. Garcia-Ojalvo, Collective stochastic coherence and synchronizability in weighted scale-free networks, New J. Phys. 16 (2014) 013036. 20 [31] M. Masoliver, N. Malik, E. Sch¨oll, A. Zakharova, Coherence resonance in a network of FitzHugh-Nagumo systems: Interplay of noise, time-delay, and topology, Chaos 27 (2017) 101102. [32] A. B. Neiman, D. F. Russell, Models of stochastic biperiodic oscillations and extended serial correlations in electroreceptors of paddlefish, Phys. Rev. E 71 (2005) 061915. [33] A. B. Neiman, D. F. Russell, Sensory coding in oscillatory electroreceptors of paddlefish, Chaos 21 (2011) 047505. [34] W. H. Nesse, L. Maler, A. Longtin, Biophysical information representation in temporally correlated spike trains, PNAS 107 (2010) 21973–21978. [35] O. Avila-Akerberg, M. J. Chacron, Nonrenewal spike train statistics: causes and functional consequences on neural coding, Exp. Brain Res. 210 (2011) 353–371. [36] T. Schwalger, B. Lindner, Patterns of interval correlations in neural oscillators with adaptation, Frontiers in Comp. Neurosci. 7 (2013) 164. [37] W. Braun, A. Longtin, Interspike interval correlations in networks of inhibitory integrate-and-fire neurons, Phys. Rev. E 99 (2019) 032402. [38] C. Bandt, B. Pompe, Permutation entropy: A natural complexity measure for time series, Phys. Rev. Lett. 88 (2002) 174102. [39] J. M. Amigo, Permutation Complexity in Dynamical Systems: Ordinal Patterns, Permutation Entropy and All That. Springer-Verlag Berlin, 2010. [40] S. P. Strong, R. Koberle, R. R. de Ruyter van Steveninck, W. Bialek, Entropy and information in neural spike trains, Phys. Rev. Lett. 80 (1998) 197–200. [41] A. Borst, F. E. Theunissen, Information theory and neural coding, Nat. Rev. Neurosci. 2 (1999) 947–957. 21 [42] A. Bahraminasab, F. Ghasemi, A. Stefanovska, P. V. E. McClintock, H. Kantz, Direction of coupling from phases of interacting oscillators: A permutation information approach, Phys. Rev. Lett. 100 (2008) 084101. [43] Z. Li, G. Ouyang, D. Li, X. Li, Characterization of the causality between spike trains with permutation conditional mutual information, Phys. Rev. E 84 (2011) 021929. [44] O. A. Rosso and C. Masoller, Detecting and quantifying stochastic and coherence resonances via information-theory complexity measurements, Phys. Rev. E 79, 040106 (2009). [45] F. Montani, R. Baravalle, L. Montangie, and O. A. Rosso, Causal information quantification of prominent dynamical features of biological neurons, Phil. Trans. Roy. Soc. A. 373, 20150109 (2015). [46] U. Parlitz, S. Berg, S. Luther, A. Schirdewan, J. Kurths, N. Wessel, Classifying cardiac biosignals using ordinal pattern statistics and symbolic dynamics, Compt. Biol. Med. 42 (2012) 319. [47] Y. Cao, W.-w. Tung, J. B. Gao, V. A. Protopopescu, L. M. Hively, Detecting dynamical changes in time series using the permutation entropy, Phys. Rev. E 70 (2004) 046217. [48] D. Arroyo, P. Chamorro, J. Amig´o, F. Rodr´ıguez, P. Varona, Event detection, multimodality and non-stationarity: Ordinal patterns, a tool to rule them all?, Eur. Phys. J. Special Topics 222 (2013) 457–472. [49] C. Quintero-Quiroz, L. Montesano, A. J. Pons, M. C. Torrent, J. Garca- Ojalvo, C. Masoller, Differentiating resting brain states using ordinal symbolic analysis, Chaos 28 (2018) 106307. [50] I. Echegoyen, V. Vera-vila, R. Sevilla-Escoboza, J. Martnez, J. Buld, Ordinal synchronization: Using ordinal patterns to capture interdependencies between time series, Chaos, Solitons & Fractals 119 (2019) 8 – 18. 22 [51] J. A. Reinoso, M. C. Torrent, C. Masoller, Emergence of spike correlations in periodically forced excitable systems, Phys. Rev. E 94 (2016) 032218. [52] R. FitzHugh, Impulses and physiological states in theoretical models of nerve membrane, Biophys. J. 1 (1961) 445. [53] J. Nagumo, S. Arimoto, S. Yoshizawa., An active pulse transmission line simulating nerve axon, Proc. IRE 50 (1962) 2061–2070. [54] M. Masoliver, C. Masoller, Sub-threshold signal encoding in coupled Fitzhugh-Nagumo neurons, Sci. Rep. 8 (2018) 8276. [55] J. A. Acebr´on, A. R. Bulsara, W.-J. Rappel, Noisy Fitzhugh-Nagumo model: From single elements to globally coupled networks, Phys. Rev. E 69 (2004) 026202. 23