1. INTRODUCTION
A single neuron reacts to the incoming stimuli in unexpected ways. As for instance, it might fires a pulse (spike) at constant rates, if the neuron operates above a threshold. The periodic firing behaviour is probably the simplest dynamical phenomena observed in the neuronal models [4], [11], [23], [33]. The variability of neuron's response to the incoming stimulus is encoded in the so-called inter-spike interval (ISI), i.e., the time-length between any two consecutive spikes [34], [35]. In a situation where the periodic regime appears, the coefficient of variations of ISIs is zero. Here, one can easily predict a sequence of spiking times once the neuron fires its first spike.
Unfortunately, the response of real neurons is far from being regular and periodic. Instead, the irregular firing activity is the most prominent feature and it might be related to brain disorders according to some experimental studies [3], [36], [37]. The irregularity is characterized by the variability on ISIs, which means that the fluctuation is not zero at all [38]. Investigating the necessary conditions in which irregular (and regular) behaviours exist is upmost important as it allows us to understand how the brain perceives the world.
Multiple sources are mainly responsible for the emergence of irregular behaviour. In the cortical neurons, for example, the intrinsic noise resulting from the fast activation of ion channels may contribute to the variability on ISIs [24]. The noise has been extensively studied for different setups and conditions. It is believed that the noise is not only slows down but also improves the reliability of neuronal firing, as well as facilitates the temporal patterning of synchronization and chaos [8], [16], [26], [27], [39], [42].
Apart from the noise effect, the irregular firing behaviour is also induced by some chemical-type currents which highlight a significant impact of synaptic transmission [7]. The real neurons communicate by sending and receiving pulse that lasts about 1-2 ms time course [11]. A recent study reports the pulse width is responsible for the irregular response of neurons in auditory systems [45]. It also has been suggested that the persistence of irregular dynamics in coupled oscillators is strongly depends on the duration of the pulse [1], [20] making it plausible for neuronal modelling.
*Corresponding author
In this paper, we investigate numerically the stochastic neural networks taking into account the pulse width as an additional parameterisation for the synaptic coupling, which has never been considered in any other similar works (see [17], [26], [27]). Moreover, compared to the previous studies, our object of study is a one-dimensional leaky integrate-and-fire (LIF) neuron [43], [44]. The model is quite simple yet the dynamics are very rich.
On top of that, the earlier works on stochastic LIF model were mostly dealing with a single neuron and concern only to the theoretical perspectives of it while less attention has been given to the emerging dynamical states in macroscopic scales [18], [30], [40], [41]. Here, We are working with a network of two noisy LIF neurons interacting via excitatory synapses to test the robustness of irregular firing activities. The synaptic transmission amongst them is modelled as a smooth function of exponential shape. Due to the randomness of stochastic input, we limit ourselves to a numerical analysis.
More precisely, in Section 2 we define the model including the Gaussian white noise and finite width pulse. In the same section, we introduce the schemes used to perform numerical simulation and the corresponding statistical quantities, namely the synchrony measure, to assess the emerging dynamical regimes. In the Section 3 we present some findings regarding the effect of noise on single neuron, which later on generalized for the two neurons interacting through exponential pulses. Finally, in Section 4 we discuss the limitations of this study and highlight some problems for the future.
2. MODEL FORMULATION
The starting point of this paper is a general stochastic system (Equation 2) introduced in [26]. Instead of electrical coupling, we incorporate the synaptic type into the model by employing the finite width pulses. The unit cell is modelled as a LIF neuron defined in [19]. We use a simple network architecture of two cells interacting via excitatory synapse with the same coupling strength as studied in [27]. More specifically, our systems (see Figure 1) consist of two identical and mutually coupled LIF neurons given by the following differential equations:
\[\dot{u} = a - u + \xi_u(t) + \frac{\mu}{2}e_u, \qquad \dot{v} = a - v + \xi_v(t) + \frac{\mu}{2}e_v,\] (1)
combined with the resetting rules: if u ≥ uth ⇒ u = ur and v ≥ vth ⇒ v = vr, respectively. The variables u and v represent the membrane potentials for each neuron, u˙ and v˙ are the corresponding derivatives with respect to time t, while a is constant current. We assume that all variables and parameters are dimensionless. The parameter µ > 0 is the coupling strength describing the intensity of excitatory postsynaptic potential. We introduce a scaling factor of 1/2 to account for the full connectivity in the present model of two interconnected neurons. In principle, the model can be generalized to N number of fully coupled oscillators, so the appropriate scaling would be 1/N [11].
Figure 1: A schematic representation for Model (1).
Following ref [1], [31], the fields eu and ev are the linear superpositions of exponential spikes emitted by
| Parameters | References | |
| Threshold: \(u_{th}, v_{th}\) | 1 | [19] |
| Reset potential: \(u_r, v_r\) | 0 | [19] |
| Constant current: a | 1.5 | [11] |
| Coupling strength: \(\mu\) | \(2 \times 10^{-3}\) | [27] |
| Noise amplitude: \(\sigma\) | [0.1, 1.4] | [26] |
| Inverse pulse width: \(\alpha\) | [20, 95] | [1] |
Table 1: The default parameters used during the simulation for Equations (1)-(2).
the sending neurons. Mathematically,
\[\dot{e}_u = -\alpha \left( e_u - \sum_k \delta \left( t - t_k^v \right) \right), \qquad \dot{e}_v = -\alpha \left( e_v - \sum_k \delta \left( t - t_k^u \right) \right), \tag{2}\]
where \(\alpha\) denotes the inverse width of the incoming excitatory pulses, and both \(t_k^u, t_k^v\) are the emission times of the kth pulse sent by the neuron u and v, respectively. Whenever u,v reaches the threshold \(u_{th}, v_{th}\), then the two neurons are reset to \(u_r, v_r\). In the limit of \(\alpha \to \infty\), the Equations (2) are simply the sum of \(\delta\) pulses, i.e.,
\[e_u = \sum_k \delta(t - t_k^v), \qquad e_v = \sum_k \delta(t - t_k^u). \tag{3}\]
The functions \(\xi_u(t)\) and \(\xi_v(t)\) are Gaussian white noise with the statistical property: \(\mathbf{E}\left[\xi(t)\right]=0\) for all time t, and the autocorrelations \(\mathbf{E}\left[\xi(t')\xi(t'')\right]=\sigma^2\delta\left(t'-t''\right)\) for \(t'\neq t''\), where \(\mathbf{E}\) denotes the expected value. Here, \(\sigma\) is the noise amplitude and \(\delta(t)\) is the delta function. For simplicity, the functions \(\xi_u(t)=\xi_v(t)=\xi(t)\), which assumes the same noisy inputs for the two neurons.
2.1. Numerical analysis
For numerical analysis, the Euler-Maruyama scheme [28] is implemented to solve simultaneously the main Equations (1)-(2). Notice that, as the integration time step \(dt \to 0\) the trajectory is smoothened since the fluctuations term vanishes with dt. For convenience, we choose \(dt = 10^{-3}\).
The initial conditions for the membrane potential are drawn randomly from a uniform distribution in a closed interval \([0,\epsilon]\) i.e., \((u(0),v(0))\in U_{[0,\epsilon]}\) while the fields are initially set to 0. Table 1 summarizes all parameters being used for the numerical simulation. Both reset potential \((u_r,v_r)\) and threshold \((u_{th},v_{th})\) have been chosen following the hypothetical values in the deterministic version of LIF neurons [19]. The spike emission of LIF neuron is achieved when the constant current is larger than the threshold, which in this case: a>1. Here we choose a=1.5 as initially set in [11]. We also assume a relatively weak coupling strength in the order of magnitude \(10^{-3}\) as introduced in [27]. Finally, the noise amplitudes and the pulse widths have been selected according to the computational experimentations [1], [26], [30].
In order to quantify the degree of synchronization in coupled systems (1)-(2), we make use of the average technique [26]
\[R = \langle \sqrt{(v-u)^2 + (e_v - e_u)^2} \rangle. \tag{4}\]
The number R is a synchronization error and the angular bracket \(\langle \star \rangle\) denotes time averaging. If the two neurons are completely synchronized then, R=0. Complete synchronization refers to a dynamical regime where all neurons within the network are firing exactly at the same time.
3. MAIN RESULTS
3.1. Single neuron
In the absence of any synaptic interaction (\(\mu = 0\)), the two neurons are indeed identical and evolve independently in time. In this case, the fields (2) become irrelevant and hence the problem of coupled
systems (1) can be reduced to a single neuron. The time evolution for a single neuron is governed by the following stochastic system:
\[\dot{y} = a - y + \xi(t), \text{ if } y \ge 1 \Rightarrow y = 0, \tag{5}\]
where y is the membrane potential. The Equation (5) is a special case of the well-known Ornstein-Uhlenbeck process. To simulate the trajectory of stochastic LIF neuron (5), it is reformulated as follows
\[dy = (a - y) dt + \sigma \sqrt{dt}x, \tag{6}\]
where x is a random number, drawn from a zero-mean Gaussian distribution with the unit variance. The same procedure of numerical integration is also adjusted for Equation (6) as in the Section 2.1.
The firing rate and coefficient of variations are statistical quantities frequently used to characterize the spiking behaviour for single neuron [11], [29]. The firing rate refers to the number of spikes divided by the length of integration time. The explicit formula is given by
\[\nu = \lim_{t \to \infty} \frac{n_t}{t},\tag{7}\]
where \(n_t\) is the number of spikes emitted by the neuron over a time interval t.
The coefficient of variations (CV) is a microscopic measure of the dynamics based on the variability of inter-spike intervals. The inter-spike interval is time difference of two consecutive spikes, i.e., \(ISI = t_{k+1} - t_k\). The CV of ISIs is estimated by a ratio
\[CV = \frac{\Delta}{\eta},\] (8)
where \(\Delta\) is standard deviation of ISIs and \(\eta\) is the corresponding mean ISIs.
When \(\sigma = 0\), meaning that no external noise is injected to the neuron, the Equation 6 is just a linear deterministic system:
\[dy = (a - y)dt. (9)\]
The only relevant parameter here is a. Note that, the neuron emits a spike if a is larger than threshold 1, otherwise there is no spike emitted. The neuron's free-noise trajectory is obtained by integrating Equation (9) with an initial condition y(0) = 0 as shown in Figure 2a. Once the membrane potential passes a threshold of 1, it is reset to 0 while at the same time, the spike is produced.
The periodic firing behaviour of the neuron implies that \(t_1\) defines a period of oscillation or ISI (for a detailed discussion on how to derive \(t_1\) explicitly see [2]). The next spikes consecutively occur at \(t_2 = 2t_1, t_3 = 3t_1, \ldots t_i = it_1\) where i is a positive integer.
Now, let us consider the dynamics of noisy neuron utilizing \(\sigma \neq 0\). The noise is characterized by a probability density function
\[f(x) = \frac{1}{\sigma\sqrt{2\pi dt}} \exp\left(-\frac{1}{2dt} \left(\frac{x}{\sigma}\right)^2\right). \tag{10}\]
Since the neuron is stimulated by the noise current with an amplitude \(\sigma=0.5\) (green colour in Figure 2c), the inter-spike interval fluctuates irregularly in time (Figure 2b). We can no longer predict the exact timing of the spikes because the spike arrivals are independent of one another. However, we may characterize the distribution of the spikes by splitting the length of integration time t into K bins such that
\[\tau = \frac{t}{K}.\tag{11}\]
If K is large (meaning that \(\tau\) becomes sufficiently small)<sup>1</sup>, the spike count in each bin is apt to be either 0 or 1. As a result, the distribution of the spikes is characterized by
\[\mathbf{E}(\{0,1\}) = \frac{n_t}{K} = \nu \tau, \qquad \mathbf{V}(\{0,1\}) = \nu \tau - (\nu \tau)^2, \tag{12}\] where \(\mathbf{E}\) and \(\mathbf{V}\) stand for the expectation and variance of the spike events, respectively. The second equality in the expected value results from the Equation 7 and 11.
<sup>1</sup>This condition implies that any two or more spikes cannot occur together, instead, either exactly there is a single spike or no spike at all.
Figure 2: (a)-(b) Time series of noise-free (\(\sigma=0\)) and noisy trajectory (\(\sigma=0.5\)) for LIF neuron (see Eqs. (5)-(6)) with a fixed a=1.5. When the membrane potential y(t) hits the threshold 1 (red dashed line), it is reset to 0. (c) The probability density for white noise when \(\sigma=0.5\) (green) and \(\sigma=1\) (black) generated from the simulation. The solid red curves correspond to f(x) in the Equation 10.
Figure 3: Firing rate \(\nu\) versus \(\sigma\), coefficient of variations CV versus \(\sigma\).
Figure 3 shows a more quantitative description of neural activity. There one can see, the firing rate \(\nu\) increases with the noise amplitude \(\sigma\) meaning that the neuron is highly sensitive to the random stimulus. Such observation is in line with the experimental study, particularly for the neurons in the cerebral cortex [6].
The same scenario is also shown when looking at the plot for the coefficient of variations CV as a function of \(\sigma\). The fluctuation of the inter-spike intervals is an effect of noisy inputs that raises by increasing the parameter \(\sigma\). Such dynamics is robust, not only in the case of LIF neuron but also found in the biophysically Hodgkin-Huxley neuronal model [21]. The noise-induced spike times variability serves as a basic mechanism for the brain functioning such as information processing and perception [9].
3.2. Noise-induced synchronization in two connected neurons
We start our analysis of the noise effect on a network of two neurons (\(\mu > 0\)) interacting through finite width pulses. All parameters are kept the same as in Section 2.1 with the initial conditions are drawn randomly in a such way that \((u(0),v(0)) \in U_{[0,1]}\). The model (1) assumes a weak coupling and so we fix the coupling strength as \(\mu=2\times 10^{-3}\). The only relevant parameters are \(\sigma\) and \(\alpha\). The simulation is performed over 10000-time units after discarding 2000-time units of transients to let enough of the systems settle onto the attractor.
Figure 4: Synchronization error R as a function of noise amplitude \(\sigma\) for different pulse widths \(\alpha^{-1}\): black triangles (\(\alpha=20\)), green circles (\(\alpha=60\)), and blue diamonds (\(\alpha=95\)). The two panels inset are the time series of the membrane potentials for \(\sigma=0.4\) and \(\sigma=1\), respectively, when \(\alpha=20\).
We turn our attention on Figure 4 where we have plotted the synchrony measure upon varying the parameter \(\sigma\) for a fixed \(\alpha\). If \(\alpha=20\), the synchronization error is approximately \(R\approx 0.5\) while the noise amplitude is within the interval range \(\sigma\in[0.1,0.59)\). In this regime, the two neurons are asynchronously firing the spike as shown by the time series of membrane potential u (solid black line) and v (red dashed line): see the upper panel inset generated for \(\sigma=0.4\).
Around \(\sigma \approx 0.6\) the synchronization error drops abruptly to 0 meaning that the two neurons are now completely synchronized. The time series analysis reveals that the ISIs fluctuate due to the presence of noise (lower panel inset). For future reference, we define \(\sigma_{\tau}\) as a transition point at which the systems collapse onto complete synchronization (CS). Formally,
\[\sigma_{\tau} = \min\{\sigma | R\left(\sigma\right) = 0\}. \tag{13}\]
It is also helpful to define a subset \(S_\tau \subset S\)
\[S_{\tau} := \{ (\alpha, \sigma_{\tau}) | R(\alpha, \sigma_{\tau}) = 0 \}, \tag{14}\] where \(S:=\{(\alpha,\sigma)\,|\,R\,(\alpha,\sigma)=0\}\). As instance, we observe the CS regime is self-sustained for \(\alpha=20\) within the range of \(\sigma\in[0.6,1.4]\), where \(\sigma_{\tau}\approx0.6\) and \(S_{\tau}=\{(20,0.6)\}\). Specifically, the time series of membrane potential show that u(t)=v(t) for each t and \(t\to\infty\): see the lower panel inset generated for \(\sigma=1\).
Upon increasing \(\alpha\) to 60 and 95, the asynchronous firing (AF) behaviour keeps persist in a finite but smaller interval range \(\sigma \in [0.1, 0.58)\) and \(\sigma \in [0.1, 0.49)\), respectively. Meanwhile, the CS emerges in a wider interval range \(\sigma \in [0.58, 1.4]\) and \(\sigma \in [0.49, 1.4]\). The Figure 5 sketches \(S_{\tau}\) for \(\alpha \in [22, 95]\) which are estimated for very long integration time. The curve \(S_{\tau}\) splits the whole plane into two regions: complete synchronization (above the curve) and asynchronous dynamics (below one).
To get a better insight into the role of noise and finite pulses to neuron dynamics, we produce a distribution of ISIs over 900000-time units for the neuron u as displayed in Figure 6. As seen, both noise and finite pulses give rise to the variations of ISIs. The width of the distributions is almost ten times wider as \(\sigma\) is increased from 0.4 to 1 for each case of \(\alpha\). Meanwhile, it is almost doubled upon shortening the pulse widths.
Figure 5: Parameter space \((\alpha, \sigma)\). The red stars correspond to \(S_{\tau}\).
Figure 6: Distribution of inter-spike intervals for \(\sigma=0.4\) (left column) and \(\sigma=1\) (right column). The vertical axis for all panels is in the logarithmic scale. Each row corresponds to different \(\alpha\): upper panel (\(\alpha=20\)), middle panel (\(\alpha=60\)) and lower panel (\(\alpha=95\)).
3.3. Coexistence regime
So far, we have explored the behaviours of the coupled neurons (1) by selecting random initial conditions \((u(0),v(0))\in U_{[0,1]}\) and found a discontinuous transition upon tuning the parameter \(\sigma\). Let us now investigate the possibility of coexistence regimes that are often accompanying the transition phenomena.
We start by fixing a small interval of width \(\epsilon\) at which the initial conditions for neuron u and v are selected whereas the fields are set to zero. Figure 7 compares the emerging scenario both for random and narrow initial conditions. In the interval range of \(\sigma \in [0.1, 0.6)\), the \(R \approx 0.5\) results from imposing random initial
Figure 7: The coexistence of complete synchronization and asynchronous firing dynamics for \(\alpha=20\). The black triangles, magenta and green crosses correspond to different initial conditions: fully random (black triangles), finite but small interval width (magenta and green crosses). The narrow initial conditions are chosen to be in the order of \(\epsilon=10^{-3}\) (magenta crosses) and \(\epsilon=10^{-2}\) (green crosses). The time series corresponds to \(\sigma=0.4\) with random initial conditions (upper panel) and narrow initial conditions (lower panel).
Figure 8: Left to right: Distribution of inter-spike intervals for neuron u when \(\alpha=20\) and \(\alpha=60\). The vertical axis is in the logarithmic scale. The colours are corresponding to different noise amplitudes: \(\sigma=0.4\) (black), \(\sigma=0.7\) (red), and \(\sigma=1\) (green). The data are all generated for the narrow initial conditions \(\epsilon=10^{-3}\).
conditions (see the black triangles). It is a clear indication of asynchronous firing dynamics. However, one can see that the complete synchronization emerges in the whole range of parameter \(\sigma\) upon choosing a small \(\epsilon = 10^{-3}\) as signalled by R = 0 (see the magenta crosses).
The scenario is further supported by the time evolution of neuron u (solid black) and v (red dashed) for \(\sigma=0.4\) in the inset panels. The upper panel correspond to the green crosses where \(\epsilon\in(10^{-3},1]\) or equivalently random initial conditions. Clearly, the two neurons fire the spike asynchronously. In comparison, the two neurons are identical for \(\epsilon=10^{-3}\) as displayed in the lower panel.
Once again, the complete synchronous state is characterized by fluctuating ISIs. We sketch the distribution of ISIs for neuron u in Figure 8. For a relatively low noise intensity (black curve), the ISIs looks homogeneous compared to the intermediate and strong (red and green curves) noises. Interestingly, decreasing the pulse width significantly adds the variability of ISIs.
4. DISCUSSION AND CONCLUSION
In this study, we have investigated the effect of noisy inputs on the firing patterns in a minimal network of two LIF neurons interacting through finite width pulses. Consistent to the earlier reports [26], [27], [14], the noise evokes the neuromodulation of synchronous and asynchronous firing dynamics. Some experimental studies show that both phenomena are robust in certain areas of the brain such as the visual cortex and hippocampus [5], [12], which might be associated with brain disease such as seizures and epilepsy [22].
We observe the existence of synchronous behaviour, as indicated by R = 0. The main important feature of the regime is the variability of ISIs, while at the same time the neurons are firing simultaneously. This study thereby suggests the dependence of the average ISI, not only on the coupling strength as highlighted in [13], but also on the noise amplitude. Given this situation, it is impossible to predict the exact timing of the spike compared to the classical case of "periodic" synchronization2 . Consequently, the computation becomes more complicated as the noise interferes with the system: one must go beyond the traditional approaches at least in the context of stability analysis [15].
As we have shown, the fluctuation of ISIs changes when the noise is varied. A sufficient strong noise causes a large excursion on the trajectory and thus greatly affects the next spiking time once the last spike has fired. This is possibly the mechanism underlying the appearance of wide ISI distributions displayed in Figure 6 (e.g., for σ = 1). Upon reducing the noise amplitudes, the two neurons desynchronize and fire the spike independently as signalled by R > 0. In the microscopic level, the regime is characterized by more uniform ISIs and hence a shorter width of ISI distributions (Figure 6, σ = 0.4), similar to that have been observed in other neuronal models [17], [25].
The novel outcome of this study is the role played by the width of exponential pulses in shaping the ISI distribution. The width of ISI distribution is lengthened as the pulse duration is shortened (see again Figure 6 and 8), which also suggesting the dependence of ISI distribution on the finite width pulses. It is unclear, however, to what extent the pulse width affects the distribution of ISIs produced by different levels of noisy inputs. This problem, in principle, can be approached mathematically by deriving an equation describing the probability distribution of ISIs by means of the Fokker-Planck equation [18], considering the pulse widths as the additional parameter.
Another observed phenomenon is the transition from complete synchronization to asynchronous dynamics. As indicated by the sharp jump of global observable variable R, the transition is discontinuous, a typical first-order transition commonly found in many living cells including neurons, as usual mechanism to respond to the external stimulus [10]. In this context, the noise should be regarded as a crucial component in the mathematical modelling for the nervous systems. The real neurons are noisy, and the noise might come from many sources such as thermal noise, ion channel propagation, synaptic transmissions, network connectivity, etc [9], [46]. Furthermore, the transition region is accompanied by a hysteresis. Upon choosing different initial conditions within a small interval range ϵ = 10−3 , we explore the emergence of bistable regimes: the new value of R = 0 (see magenta crosses in Figure 7) indicates a complete synchronous behaviour.
Our findings highlight the key points for future reference. First, even in simple stochastic LIF neural networks, the dynamical regimes are intriguing and are in line with the studies in a higher dimensional system, such as Hindmarsh-Rose and Morris-Lecar neurons (see again [26], [27]). We hypothesize the overall scenario emerging in this model is likely to persist in the larger system sizes. Secondly, the finite pulses is matter: shortening the pulse width adds the variability of ISIs in stochastic neurons. If confirmed analytically, these statistics could lead to well-established computational models and thus better understanding of irregular activity in the brain.
Finally, it is important to note some limitations of this study. We are working with a minimal network of two globally coupled identical neurons, while the brain consists of billions of neurons having more complex structures and intrinsic properties. The works herein can be extended to large-scale networks in the assumption of two interacting populations of excitatory and inhibitory neurons, which will be our next object of study. This study is also limited to test numerically the impacts of noise and finite pulses on neuron dynamics (see again Figure 5). We argue that the asynchronous firing dynamics vanishes for a large α, therefore further studies are required.
2 "Periodic" synchronization refers to the regimes where the fluctuation of ISIs of the whole networks is zero on average [1], [32].
ACKNOWLEDGEMENT
This work was partially supported by the Indonesian Mathematical Society (IndoMS) through a visiting research program, grant No. 028/Pres/IndoMS/SP/VIII/2022.