Sub-threshold signal encoding in coupled FitzHugh-Nagumo neurons

Despite intensive research, the mechanisms underlying the neural code remain poorly understood. Recent work has focused on the response of a single neuron to a weak, sub-threshold periodic signal. By simulating the stochastic FitzHugh-Nagumo (FHN) model and then using a symbolic method to analyze the firing activity, preferred and infrequent spike patterns (defined by the relative timing of the spikes) were detected, whose probabilities encode information about the signal. As not individual neurons but neuronal populations are responsible for sensory coding and information transfer, a relevant question is how a second neuron, which does not perceive the signal, affects the detection and the encoding of the signal, done by the first neuron. Through simulations of two stochastic FHN neurons we show that the encoding of a sub-threshold signal in symbolic spike patterns is a plausible mechanism. The neuron that perceives the signal fires a spike train that, despite having an almost random temporal structure, has preferred and infrequent patterns which carry information about the signal. Our findings could be relevant for sensory systems composed by two noisy neurons, when only one detects a weak external input.


Results
We simulate the coupled FHN neurons as described in Methods, with a periodic sub-threshold signal that is applied to one of the neurons, referred to as neuron 1. Figure 1 displays the voltage-like variable of neuron 1, u 1 , in different situations. When there is no noise, no signal and no coupling, the neuron is in the rest state and when the sub-threshold signal is applied, u 1 displays small sub-threshold oscillations [panel (a)]; when noise is added, noise-induced spikes are observed, which carry information about the applied sub-threshold signal [panel (b)]; and when the coupling to neuron 2 is added, a noticeable effect is the increase of the firing rate [panel (c)]. The differences that are qualitatively observed in these time series are going to be quantitatively addressed by using the methods of analysis presented in Methods.
As we are interested in the encoding of weak signals, we first have to distinguish between a sub-threshold and a super-threshold signal. The first one refers to a signal which, in the absence of noise, it does not induce any spike [u 1 displays small oscillations, as in Fig. 1(a)], while the second one is a signal that is strong enough to induce spikes. A periodic signal can be either sub-threshold or super-threshold depending on both, the period and the amplitude. Thus, to identify the parameters where the signal is sub-threshold, in Fig. 2 we plot in color code the spike rate (i.e., the inverse of the mean ISI, 1/〈I〉), as a function of a 0 and T. In panel (a) neuron 1 is isolated (σ 2 = 0), while in panel (b) it is coupled to neuron 2 (σ 1 = σ 2 = 0.05).
When the neuron is uncoupled, for large amplitude and/or small period the signal is super-threshold, otherwise is sub-threshold. When the neuron is coupled to neuron 2, we note that the super-threshold region is slightly larger in the parameter space (a 0 , T), as compared to the uncoupled case.
When we include noise, Fig. 2(c) and (d), we first note that in the super-threshold region (yellow) the spike rate does not change significantly (it is about the same as for D = 0). This is due to the fact that in this region the spikes are mainly induced by the signal.
In contrast, in the sub-threshold region, comparing the uncoupled (panel c) and the coupled (panel d) situations, we note that coupling significantly increases the spike rate (it almost doubles). Therefore, in this region coupling plays the role of an extra source of noise.
Having identified the sub-threshold region in the parameter space (a 0 , T), we next turn our attention to the influence of the coupling coefficients. Figure 3 displays the spike rate as a function of σ 1 and σ 2 in different situations. In panel (a) there is no signal and no noise. We observe that when both |σ 1 | and |σ 2 | are large enough, the coupling induces spikes. Positive coupling coefficients result in a higher spike rate, in comparison with negative coefficients. In panel (b), the noise is still zero but a weak signal is applied. Because the signal is subthreshold [a 0 = 0.05 and T = 10, which are in the sub-threshold region in Fig. 2(a) and (b)], we note only small variations with respect to panel (a).
In Fig. 3(c) and (d) noise is included; in (c) there is no signal while in (d) the weak signal is applied. To show how the spike rate changes with the coupling, Fig. 3(c) and (d) display the relative variation of the spike rate (with respect to the spike rate when neuron 1 is uncoupled). Without signal (panel c), positive coupling coefficients result in larger spike rate as compared to negative ones, however, when the signal is applied (panel d) these differences are washed out. The vertical line in panels (c) and (d) is due to the fact that when σ 1 = 0 neuron 1 is uncoupled from neuron 2, and thus its spike rate does not depend of σ 2 .
In order to limit the number of parameters, in the following we assume σ 1 = σ 2 = σ and fix (unless otherwise stated) σ = 0.05, a 0 = 0.05 and T = 10. For these parameters the signal and the coupling act as sub-threshold perturbations: without noise neuron 1 does not fire spikes.
To further characterize the role of noise, Fig. 4 displays the mean ISI, 〈I〉, as a function of noise intensity for different periods of the applied signal (in the Supplementary Information we analyze the shape of the ISI distribution). In panel (a) σ = 0, while in panel (b), σ = 0.05. For both cases there is clearly a noise dominated regime, where 〈I〉 is the same, regardless of the coupling and of the period of the signal. In contrast, for low noise levels  the coupling and the period affect the 〈I〉. In panel (a) (σ = 0) we can also compare the mean ISI when the signal is applied (solid symbols indicate a 0 ≠ 0 and different periods) and when the signal is not applied (empty circles): we see that, when a 0 ≠ 0 the neuron fires at lower noise intensities as compared to a 0 = 0. Comparing panel (a) with panel (b) (σ = 0.05) we note that when neuron 1 is coupled to neuron 2, it starts firing at even lower noise intensities.
Regarding the role of the period of signal, when the noise level is low, the larger T is, the larger 〈I〉 is. There is a linear relation, as shown in Fig. 4(c) and (d), which holds for both, the coupled and the uncoupled cases. For stronger noise, 〈I〉 remains constant when increasing T.
Noise-induced regularity in the spike train [30][31][32] is characterized in panels (e) and (f), where the normalized standard deviation of the ISI distribution, R, is plotted against the noise intensity for different T, without and with coupling, respectively. In both panels, two minimums are observed. Whereas the first one indicates stochastic resonance [33][34][35] , as it occurs when ~〈 〉 T I , the second one reveals the coherence resonance phenomenon 30,36 , which is independent from the period of the signal. It occurs for an intermediate value of the noise amplitude for which noise-induced oscillations become most coherent. For some periods T a maximum appears for very small values of the intensity of the noise. Such maxima are a signature of anticoherence resonance 37 .
After having characterized the effects of the weak signal, of the coupling, and of the noise in the neuron's spike rate and in the regularity of the spikes, we next turn our attention to the timing of the spikes. We use symbolic ordinal analysis (see Methods) to investigate the possible presence, in the ISI sequence, of over expressed and of less expressed spike patterns.
We begin by considering the situation in which no signal is applied and analyze the effect of increasing the noise level or the coupling strength: Fig  respectively. We note that neither the noise nor the coupling induce temporal order in the spike sequence (as all the probabilities are within the blue region that indicates values consistent with equal probabilities). When the signal is applied, panels (c) and (d), we note that there is temporal order in the spike sequence, as the ordinal probabilities reveal the presence of over expressed and less expressed spike patterns (the probabilities are not in the blue region and thus, are not consistent with the uniform distribution). Moreover, we note that the variation of the probabilities with D or σ is qualitatively similar.
Next, we analyze how the coupling coefficients affect the encoding of the signal features (the amplitude and period): we compare how the ordinal probabilities vary with a 0 and T, when neuron 1 is isolated [ Fig. 6(a) and (c)] and when it is coupled to neuron 2 [ Fig. 6(b) and (d)]. In both cases, when a 0 increases (within the sub-threshold region) the probabilities monotonically increase or decrease. This variation is consistent with the results reported in 22 . While in 22 the sub-threshold signal was applied to the slow variable, v, here it is applied to the fast variable, u. In both cases, the probabilities encode information of the amplitude of the signal. Nevertheless, coupling to neuron 2 changes the preferred and infrequent patters, i.e., modifies the temporal order in the spike sequence. For instance, for σ = 0.05 the probability of the ordinal pattern 012 monotonically increases with a 0 , whereas for σ = 0.05 monotonically decreases. In panels (b) and (d) we note that, with or without coupling, the preferred and infrequent patterns depend on the period of the signal, confirming the results reported in 22 .
Comparing Fig. 6(c) and (d) we note that the coupling can either improve or degrade the signal encoding with respect to the uncoupled situation: for T = 6 and T = 10 with coupling (panel d) the probabilities have extreme values (maximum or minimum depending of the ordinal pattern), and thus, the coupling enhances the signal encoding. In contrast, for T ~ 17 with coupling (panel d) all the probabilities are close to the blue region (while without coupling they are not), which means that with coupling the probabilities do not encode information of the signal.
Next we investigate if there is an optimal combination of the signal period, T, and the coupling coefficients, σ 1 and σ 2 , for signal encoding. To quantify the information content of the spike train (represented by symbolic ordinal patterns) we calculate the entropy computed from the probabilities of the ordinal patterns (known as permutation entropy, = −∑ H p p log i i i 25 ) and normalize the entropy to its maximum value, H max = −logL! with L! being the possible number of patterns (see Methods). Figure 7(a)-(c) display the normalized permutation entropy in color code as a function of σ 1 and σ 2 for T = 6, T = 10 and T = 14, respectively. We observe values very close to 1, which indicate that the timing of the spikes is almost random (the ordinal probabilities are almost equal). This is expected as the signal parameters and the coupling strengths are sub-threshold, i.e., the spiking activity is due to the presence of noise (without noise, the neuron displays sub-threshold oscillations). However, for T = 10 (panel b) we see that when σ 1 σ 2 > 0 the entropy slightly decreases, which indicates that there are more and less expressed patterns in the spike sequence, i.e., the spike sequence carries information about the signal.
It is interesting to compare the results obtained with nonlinear ordinal analysis, with those obtained with linear analysis. Linear correlations between inter-spike intervals are detected by the serial correlation coefficients (SCCs, see Methods). In Fig. 8 the ordinal probabilities and the SCCs are plotted vs. the mean ISI, 〈I〉, which is tuned by changing the noise strength [increasing D decreases 〈I〉 as shown in Fig. 4(a) and (b)]. We see that when the noise is strong (i.e. small 〈I〉), the ordinal probabilities are outside the blue region and thus capture temporal ordering in the ISI sequence; in contrast, C 1 and C 2 (that are close to zero) do not capture linear correlations.  Another relevant issue to discuss is how the coupling terms are implemented. While we have presented model simulations where the terms σ 2 u 1 and σ 1 u 2 couple neuron 1 to neuron 2 and vice-versa 38 (see Methods), we have also simulated the model with (i) the coupling in the recovery-like variable (i.e., σ 2 v 1 and σ 1 v 2 added to the rate equations of v 2 and v 1 respectively) and (ii) with differential coupling (i.e., σ(u 1 − u 2 ) and σ(u 2 − u 1 ) added to the rate equations of u 1 and u 2 respectively). We have consistently found that the probabilities of the ordinal patterns vary with both, the period and the amplitude of the signal, in a similar way as with with non diffusive coupling (see Fig. 9).

Discussion
We have studied two coupled neurons with a weak sub-threshold periodic signal applied to one of them. We have analyzed how the firing activity of the neuron that perceives the signal encodes the signal information, and the role of another neuron that does not perceive the weak signal. We have shown that when the neuron that perceives the signal is coupled to the second neuron, the spike rate increases and the noise level needed for firing spikes decreases, with respect to the uncoupled neuron. We have used symbolic ordinal analysis to investigate temporal ordering in the timing of the spikes fired by the neuron that perceives the signal. We have shown that the spike sequence has over expressed and less expressed ordinal patterns whose probabilities carry information about the features of the signal (the amplitude and the period). We have also shown that, when the noise is strong, the ordinal probabilities can still encode information about the weak signal, which is not encoded in the spike rate (that is independent of the period of the signal) and is not detected by linear correlation analysis (as the serial correlation coefficients at lags 1 and 2 vanish).
Clearly, it is crucial that the neuron that perceives the signal not only encodes the information, but also, transfers the information. In order to investigate information transfer, we plan to analyze, in future work, how the spike sequence of the second neuron (that does not perceive the signal) encodes the information of the signal perceived by the first neuron. In this sense, it is important to compare the dynamical behavior of the second neuron, with  and without the applied weak signal, in order to determine which specific property of the spike train (spike rate, SCCs, ordinal probabilities) carry information about the features of the weak signal (amplitude and period) that is perceived by the first neuron. It would also be interesting to analyze the encoding of more complicated signals, for example, a weak signal with two different frequencies, and compare with the phenomenon of vibrational resonance 39 .
Our findings could be relevant for neuronal sensory systems composed by coupled noisy neurons, when only one is affected by external inputs. The encoding mechanism demonstrated here, by which the period and the amplitude of the applied sub-threshold signal are encoded in the values of the ordinal probabilities, is very slow if the probabilities are computed from the spike train of a single neuron, because a large number of spikes are needed in order to compute the patterns' probabilities. However, if the encoding is performed by neuronal populations, then, the probabilities can be computed from the spikes of many neurons, and in this case, only few spikes per neuron would be enough to compute the probabilities. This ensemble-based encoding mechanism would allow fast encoding of time-varying signals. Ongoing work is devoted to understand the robustness of the proposed signal encoding mechanism when the neuron (or neurons) that perceive the signal is (are) coupled in a small modular network. We are also studying the role of non-Gaussian noise (Poisson), and of synaptic coupling (excitatory or inhibitory). If the encoding mechanism is indeed robust, we then plan to investigate information transmission in large neuronal populations [40][41][42][43][44] .

Methods
Model. We consider two identical FitzHugh-Nagumo neurons 23,24 (in the Supplementary Information we present simulations of non-identical neurons), mutually coupled as in 38 , with a periodic signal applied to one of them (referred to as neuron 1):     The dimensionless variables u i and v i are a fast variable that represents the voltage of the membrane, and a recovery-like variable that represents the refractory properties of the membrane (slow variable); a and ε are parameters that control the spiking activity of the uncoupled neurons. The coupling terms σ 2 u 1 and σ 1 u 2 mimic synaptic currents from neuron 1 to neuron 2 and vice-versa 38 . The signal has amplitude a 0 and period T. The noise is modeled with statistically independent Gaussian white noise terms [〈ξ i (t)ξ i (t′)〉 = δ(t − t′) and 〈ξ i (t)ξ− j (t)〉 = δ(i − j)] and the noise level, D, is the same for both neurons.
The values of the parameters, a = 1.05 and ε = 0.01, are chosen such that, when D = 0 and σ 1 = σ 2 = 0, the neurons are in the excitable regime: each neuron resides in a stable state (rest state) unless it is perturbed. If a strong enough perturbation occurs, the neuron leaves the rest state and after firing a spike, it returns to the rest state. Then, a refractory period follows during which another perturbation will not trigger a spike.
The equations are integrated, starting from random initial conditions, using the Euler-Maruyama method with an integration step of dt = 10 −3 . The signal parameters, a 0 and T, and the coupling coefficients, σ 1 and σ 2 , are varied within the "sub-threshold" region of the parameter space: without noise the voltage-like variables u 1 and u 2 display only small oscillations [see Fig. 1(a)]. For each set of parameters, the voltage-like variable of the neuron that receives the signal, u 1 , is analyzed and the ISI sequence is computed, {I i ; I i = t i+1 − t i } with t i defined by the condition u 1 (t i ) = 0 considering only the ascensions.
To compute the mean ISI and the coefficient R (see below) time series with a minimum number of 100 spikes are generated (as this is sufficient to estimate the mean values of the ISI distribution), while to compute the ordinal probabilities, time series with at least 10000 spikes are generated. This is because a large number of ordinal patterns are needed in order to determine if their probabilities are consistent or not with the uniform distribution 29 .
Analysis of ISI sequences. The regularity of the ISI sequence is characterized by the coefficient R 30 : where 〈I〉 is the mean value of the ISI distribution.
Correlations between ISIs are characterized by the serial correlation coefficients (SCCs): where j is an integer number. SCCs are a standard tool to analyze spike trains, however, they only capture linear correlations. In contrast, a symbolic methodology known as ordinal analysis 25 has been demonstrated to be well suited for detecting nonlinear correlations in spike trains 22,26,29  spike train of neuron 1: the ISI sequence {I 1 , ..., I i , ..., I N } is transformed into a sequence of ordinal patterns, which are defined by the relative order of L consecutive ISI values.
Once the length L of the ordinal patterns is defined, for each interval I i the subsequent L − 1 intervals are considered and compared. The total number of possible order relations (i.e., ordinal patterns of length L) is then equal to the number of permutations L!. If we set L = 2 we have only two patterns: 12 and 21 for I 1 < I 2 and I 1 > I 2 , respectively; if we set L = 3, we have 3! = 6 possible ordinal patterns.
The symbolic sequence of ordinal patterns is computed using the function perm indices defined in 45 . Then, the ordinal probabilities are estimated as p i = N i /M where N i denotes the number of times the i-th pattern occurs in the sequence, and M denotes the total number of patterns. If the patterns are equi-probable one can infer that there are no preferred order relations in the timing of the spikes. On the other hand, the presence of frequent (or infrequent) patterns will result into a non-uniform distribution of the ordinal patterns. A binomial test will be used to analyze the significance of preferred and infrequent patterns: if all the ordinal probabilities are within the interval [p − 3σ p , p + 3σ p ] (with p = 1/L! and σ = − p p M (1 )/ p ), the probabilities are consistent with the uniform distribution, else, there are significant deviations which reveal the presence of over expressed and less expressed patterns.
Here we use L = 3, which allows to investigate order relations among three ISI (i.e., four consecutive spike times). This choice is motivated by the fact that the signal is sub-threshold, i.e., the firing activity of neuron 1 is driven by noise (without noise, there are no spikes). Therefore, only short ISI correlations are expected in the spike train.