The stochastic Fitzhugh-Nagumo neuron model in the excitable regime embeds a leaky integrate-and-fire model.

In this paper, we provide a complete mathematical construction for a stochastic leaky-integrate-and-fire model (LIF) mimicking the interspike interval (ISI) statistics of a stochastic FitzHugh-Nagumo neuron model (FHN) in the excitable regime, where the unique fixed point is stable. Under specific types of noises, we prove that there exists a global random attractor for the stochastic FHN system. The linearization method is then applied to estimate the firing time and to derive the associated radial equation representing a LIF equation. This result confirms the previous prediction in Ditlevsen and Greenwood (J Math Biol 67(2):239-259, 2013) for the Morris-Lecar neuron model in the bistability regime consisting of a stable fixed point and a stable limit cycle.

Introduction

Mathematical modeling has emerged as an important tool to handle the overwhelming structural complexity of neuronal processes and to gain a better understanding of their functioning from the dynamics of their model equations. However, the mathematical analysis of biophysically realistic neuron models such as the 4-dimensional Hodgkin–Huxley (HH) (1952) and the 2-dimensional Morris–Lecar (ML) (1981) equations is difficult, as a result of a large parameter space, strong nonlinearities, and a high dimensional phase space of the model equations. The search for simpler, mathematically tractable (small parameter space, weaker nonlinearities, low dimensional phase space) neuron models that still capture all, or at least some important dynamical behaviors of biophysical neurons (HH and ML) has been an active area of research.

The efforts in this area of research have resulted in easily computable neuron models which mimic some of the dynamics of biophysical neuron models. One of the resulting models is the 2-dimensional FitzHugh–Nagumo (FHN) neuron model (FitzHugh 1961). The FHN model has been so successful, because it is at the same time mathematically simple and produces a rich dynamical behavior that makes it a model system in many regards, as it reproduces the main dynamical features of the HH model. In fact, the HH model has two types of variables, and each type then is combined into a single variable in FHN: The (V, m) variables of HH correspond to the v variable in FHN, whose fast dynamics represents excitability; the (h, n) variables correspond to the w variable, whose slow dynamics represents accommodation and refractoriness.

The fact that the FHN model is low dimensional makes it possible to visualize the solution and to explain in geometric terms important phenomena related to the excitability and action potential generation mechanisms observed in biological neurons. Of course, this comes at the expense of numerical agreement with the biophysical neuron models (Yamakou 2018). The purpose of the model is not a close match with biophysically realistic high dimensional models, but rather a mathematical explanation of the essential dynamical mechanism behind the firing of a neuron. Moreover, the analysis of such simpler neuron models may lead to the discovery of new phenomena, for which we may then search in the biological neuron models and also in experimental preparations.

There is, however, an even simpler model than FHN, the leaky integrate-and-fire model (LIF). This is the simplest reasonable neuron model. It only requires a few basic facts about nerve cells: they have membranes, they are semipermeable, and they are polarizable. This suffices to deduce a circuit equivalent to that of the membrane potential of the neuron: a resistor-capacitor circuit. Such circuits charge up slowly when presented with a current, cross a threshold voltage (a spike), then slowly discharge. This behavior is modeled by a simple 1D equation together with a reset mechanism: the leaky integrate-and-fire neuron model equation (Gerstner and Kistler 2002). Combining sub-threshold dynamics with firing rules has led to a variety of 1D leaky integrate-and-fire descriptions of a neuron with a fixed membrane potential firing threshold (Gerstner and Kistler 2002; Lansky and Ditlevsen 2008) or with a firing rate depending more sensitively on the membrane potential (Pfister et al. 2006). In contrast to dimensional neuron models, , such as the HH, ML, and FHN models, the LIF class of neuron models is less expensive in numerical simulations, which is an essential advantage when a large network of coupled neurons is considered.

Noise is ubiquitous in neural systems and it may arise from many different sources. One source may come from synaptic noise, that is, the quasi-random release of neurotransmitters by synapses or random synaptic input from other neurons. As a consequence of synaptic coupling, real neurons operate in the presence of synaptic noise. Therefore, most works in computational neuroscience address modifications in neural activity arising from synaptic noise. Its significance can however be judged only if its consequences can be separated from the internal noise, generated by the operations of ionic channels (Calvin and Stevens 1967). The latter is channel noise, that is, the random switching of ion channels. In many papers channel noise is assumed to be minimal, because typically a large number of ion channels is involved and fluctuations should average out, and therefore, the effects of synaptic noise should dominate. Consequently, channel noise is frequently ignored in the mathematical modeling. However, the presence of channel noise can also greatly modify the behavior of neurons (White et al. 2000). Therefore, in this paper, we study the effect of channel noise. Specifically, we add a noise term to the right-hand side of the gating equations (the equation for the ionic current variable).

In the stochastic model, the deterministic fixed point is no longer a solution of the system. The fixed point necessarily needs to vary and adapt to the noise. To account for this, in the theory of random dynamical systems, the notion of a random dynamical attractor was developed as a substitute for deterministic attractors in the presence of noise. In the first part of this paper, we therefore prove that our system admits a global random attractor, for both additive and multiplicative channel noises. This can be seen as a theoretical grounding of our setting.

In Ditlevsen and Greenwood (2013), it was shown that a stochastic LIF model constructed with a radial Ornstein–Uhlenbeck process is embedded in the ML model (in a bistable regime consisting of a fixed point and limit cycle) as an integral part of it, closely approximating the sub-threshold fluctuations of the ML dynamics. This result suggests that the firing pattern of a stochastic ML can be recreated using the embedded LIF together with a ML stochastic firing mechanism. The LIF model embedded in the ML model captures sub-threshold dynamics of a combination of the membrane potential and ion channels. Therefore, results that can be readily obtained for LIF models can also yield insight about ML models. In the second part of this paper, we here address the problem to obtain a stochastic LIF model mimicking the interspike interval (ISI) statistics of the stochastic FHN model in the excitable regime, where the unique fixed point is stable. Theoretically, we obtain such a LIF model by reducing the 2D FHN model to the one dimensional system that models the distance of the solution to the random attractor as shown in the first part of the paper. In fact, we show that this distance can be approximated to the fixed point, up to a rescaling, as the Euclidean norm of the solution of the linearization of the stochastic FHN equation along the deterministic equilibrium point, and hence the LIF model is approximated by the equation for . An action potential (a spike) is produced when exceeds a certain firing threshold . After firing the process is reset and time is back to zero. The ISI is identified with the first-passage time of the threshold, , which then acts as an upper bound of the spiking time of the original system. By defining the firing as a series of first-passage times, the 1D radial process together with a simple firing mechanism based on the detailed FHN model (in the excitable regime), the firing statistics is shown to reproduce the 2D FHN ISI distribution. We also show that and share the same distribution.

The rest of the paper is organized as follows: Sect. 2 introduces the deterministic version of the FHN neuron model, where we determine the parameter values for which the model is in the excitable regime. In Sect. 3, we prove the existence of a global random attractor of the random dynamical system generated by the stochastic FHN equation; and furthermore derive a rough estimate for the firing time using the linearization method. The corresponding stochastic LIF equation is then derived in Sect. 4 and its distribution of interspike-intervals is found to numerically match the stochastic FHN model.

The deterministic model and the excitable regime

In the fast time scale t, the deterministic FHN neuron model iswhere is the activity of the membrane potential and is the recovery current that restores the resting state of the model. I is a constant bias current which can be considered as the effective external input current. is a small singular perturbation parameter which determines the time scale separation between the fast t and the slow time scale . Thus, the dynamics of is much faster than that of . and are parameters.

The deterministic critical manifold defining the set of equilibria of the layer problem associated to Eq. (2.1) (i.e., the equation obtained from Eq. (2.1) in the singular limit , see Kuehn (2015) for a comprehensive introduction to slow-fast analysis), is obtained by solving for w. Thus, it is given byWe note that for Eq. (2.1), coincides with the v-nullcline (the red curve in Fig. (1)). The stability of points on as steady states of the layer problem associated to Eq. (2.1) is determined by the Jacobian scalar . This shows that on the critical manifold, points with are stable while points with are unstable. It follows that the branch is stable, is unstable, and is stable.

The set of fixed points which define the resting states of the neuron is given byThe sign of the discriminant , determines the number of fixed points. can therefore intersect the w-nullcline () at one, two or three different fixed points. We assume in this paper that , in which case we have a unique fixed point given bywhereHere, we want to consider the neuron in the excitable regime (Ditlevsen and Greenwood 2013). A neuron is in the excitable regime when starting in the basin of attraction of a unique stable fixed point, an external pulse will result into at most one large excursion (spike) into the phase space after which the phase trajectory returns back to this fixed point and stays there (Izhikevich 2007).

In order to have Eq. (2.1) in the excitable regime, we choose and such that (i.e., a unique fixed point) and such that the Jacobian (the linearization matrix M) of Eq. (2.1) at the fixed point has a pair of complex conjugate eigenvalueswith negative real part (i.e., a stable fixed point). In that case, is the only stationary state and there is no limit cycle of system (2.1). In other words, is the global attractor of the system (Izhikevich 2007). Moreover, to apply the averaging technique (Baxendale and Greenwood 2011), it is necessary that , we therefore use throughout this paper the following parameters of the system: so that is the unique stable fixed point and . Figure (1) shows the neuron in the excitable regime. Notice that although every trajectory finally converges to the fixed point, only a small change in the location of the starting point will result in different behavior of the trajectories (see the blue and purple curves).

The stochastic model

We consider this stochastic FHN modelwhere the deterministic fields f and g are given in Eq. (2.1). There are two important cases: either (additive channel noise) or (multiplicative channel noise). stands for the Stratonovich stochastic integral with respect to the Brownian motion .

Figure 2 shows the phase portraits of Eq. (3.1) starting with the initial condition , which is in the vicinity of the stable fixed point. Given an initial condition close to the stable fixed point , the trajectory of the stochastic system might first rotate around the stable fixed point but then the noise may trigger a spike, that is, a large excursion into the phase space, before returning to the neighbourhood of the fixed point; the process repeats itself leading to alternations of small and large oscillations. A similar behavior can be observed when the deterministic system with an additional limit cycle is perturbed by noise (as seen in the bistable system Ditlevsen and Greenwood 2013).

Fig. 2

Random trajectory of Eq. (3.1) in the excitable regime with chosen parameters of the system and the initial condition for both additive and multiplicative noise (we use the StochasticRungeKutta method in Mathematica with the real time and the step size )

Figure 3 shows that the spiking frequency increases as the amplitude of the noise increases. For a fixed simulation time , the system spikes only rarely, if at all, when the amplitude , but spikes more frequently when increases. This is similar for multiplicative noise.

Fig. 3

The components (left column: , right column: ) of a random trajectory of Eq. (3.1) in the excitable regime with chosen parameters of the system and the initial condition for additive noise with

Let and be the drift and diffusion coefficients of (3.1). The stochastic system is then of the formwhere for additive noise and for multiplicative noise. It is easy to check that F is dissipative in the weak sense, i.e.whereOn the other hand, we havefor multiplicative noise, while for additive noise, so H is globally Lipschitz continuous.

The existence of a random attractor

In the sequel, we are going to prove that there exists a unique solution of (3.1) and the solution then generates a so-called random dynamical system (see e.g. Arnold 1998, Chapters 1–2).

More precisely, let be a probability space on which our Brownian motion is defined. In our setting, can be chosen as , the space of continuous real functions on which are zero at zero, equipped with the compact open topology given by the uniform convergence on compact intervals in , as , the associated Borel--algebra and as the Wiener measure. The Brownian motion can then be constructed as the canonical version .

On this probability space we construct a dynamical system as the Wiener shiftThen satisfies the group property, i.e. for all , and is -preserving, i.e. for every , . The quadruple is called a metric dynamical system.

Given such a probabilistic setting, Theorem 3.1 below proves that the solution mapping defined by is a random dynamical system satisfying and the cocycle propertyTo investigate the asymptotic behavior of the system under the influence of noise, we shall first check the effect of the noise amplitude on firing. Under the stochastic scenario, the fixed point is no longer the stationary state of the stochastic system (3.1). Instead, we need to find the global asymptotic state as a compact random set depending measurably on such that A is invariant under , i.e. , and attracts all other compact random sets in the pullback sense, i.e.where d(B|A) is the Hausdorff semi-distance. Such a structure is called a random attractor (see e.g. Crauel et al. 1997 or Arnold 1998, Chapter 9).

The following theorem ensures that the stochastic system (3.1) has a global random pullback attractor. The proof is provided in the “Appendix”.

Theorem 3.1

There exists a unique solution of (3.2) which generates a random dynamical system. Moreover, the system possesses a global random pullback attractor.

Theorem 3.1 shows that every trajectory would in the long run converge to the global random attractor. The structure and the inside dynamics of the global random attractor are still open issues which might help understand the firing mechanism.

The normal form at the equilibrium point

One way to study the dynamics of the stochastic system (3.1) is through its linearization. Therefore, in this section, we shall study the dynamics of (3.1) in a small vicinity of the fixed point . To do that, consider the shift system w.r.t. the fixed point which has the formwith initial point , where is the linearized matrix of F at , is the nonlinear term such thatfor an increasing function , , which implies that . Since is either a constant or a linear function, we prove below that system (3.8) can be well approximated by its linearized system

Theorem 3.2

Given and equations (3.8), (3.9), define the stopping time . Then there exists a constant C independent of r such that for any , the following estimates hold

For additive noise

For multiplicative noise

The proof is provided in the “Appendix”. In practice we can even approximate (3.8) by the following linear system with additive noiseBy the same arguments as in the proof of Theorem 3.2, we can prove the following estimatefor the same stopping time .

Another comparison between the processes and can be obtained by using power spectral density estimation (see, for example, Fan and Yao 2003, Chapter 7). In Fig. 4, the estimated spectral densities of the shifted original and the linearized process are plotted. The spectral densities are estimated from paths started from 0 to 50 ms of subthreshold fluctuations, and scaled to have the same maximum at 40.

Fig. 4

The process (3.8) and the process  (3.9) with the chosen parameters of the system, and the same starting point are compared by using the power spectral density. Their spectrum densities are well approximated

The embedded LIF model

In this section, we present two constructive methods to obtain 1-D LIF models corresponding to the stochastic FHN in the excitable regime in Eq. (3.1). The first method follows Baxendale and Greenwood (2011) (see also Ditlevsen and Greenwood 2013) by constructing the so-called radial Ornstein–Uhlenbeck equation. More precisely, we rewrite the linearized system (3.9) in the formwhere and is a 2-D standard Brownian motion. For chosen parameters, M has a pair of complex conjugate eigenvalues with . By transformation with we obtainwhereWe note that , therefore, by applying the technique of time average from (Baxendale and Greenwood 2011, Theorem 1), can be approximated by an Ornstein-Uhlenbeck process up to a rotation, i.e.where , the rotationand is the unique solution of the 2-D SDEwith the initial value . Therefore, can be approximated by which by Ito calculus satisfies the SDEThe second method is to consider in polar coordinates withwhere . Its norm and its angle satisfyBy the averaging technique from (Baxendale and Greenwood 2011, Theorem 1) one can approximate , henceThus, by using the averaging technique, we proved that both Eqs. (4.3) and (4.4) are good approximations of the radial process . This can also be tested by using the power spectral density estimation (see Fig. 5).

Fig. 5

The norm of the process (4.2), the process  (4.3), and the process  (4.4) with the chosen parameters of the system, and the same starting point are compared by using the power spectral density. Their spectrum densities are well approximated (color figure online)

Firing mechanism

A spike in Eq. (3.1) occurs when there is a transition of a random trajectory from the vicinity of the stable fixed point located on the left stable part of to its right stable part and back to the vicinity of . This spike happens almost surely when a random trajectory with the starting point in the vicinity of crosses the threshold line . From the phase space of Eq. (3.1) (see Fig. 2), the probability of a spike increases as the starting point moves farther away from .

In order to construct the firing mechanism of Eq. (4.3) matching that of Eq. (3.1), we will calculate the conditional probability that Eq. (3.1) fires given that the trajectory crosses the line . Denote by with for , then the distance between the equilibrium and is . The value can be considered as the distance between the fixed point and the separatrix (see also Fig. 1) along L. For a given pair (), a short trajectory starting in was simulated from (3.1), it was recorded whether a spike occurred (crossing the threshold ) in the first cycle of the stochastic path around . This was repeated 1000 times and we counted the ratio of the number of spikes, denoted by , which is an estimate for the conditional probability of firing . The estimation was, furthermore, repeated for .

Fig. 1

The critical manifold (red curve) and the w-nullcline (green line) intersect at the unique and stable fixed point . Two deterministic trajectories are shown, the purple curve starts at and the blue curve starts at . Parameters of the system and the real time for trajectories (color figure online)

From the numerical simulation, for each , the estimate of the conditional probability is close to zero when we start in the immediate neighborhood of the stable fixed point and close to one when we start at the , i.e., sufficiently far from the fixed point. Theses estimates appear to depend in a sigmoidal way on the distance from the stable fixed point. Therefore we assumed the conditional probability of firing to be of the formThe parameters a and b then are estimated by using a non-linear regression from the above simulation data and are plotted in Fig. 6 for some different values of the noise amplitude , and 0.015. We see that the family of estimates, , fits the fitted curve quite well for each value of . Regression estimates are reported in Table 1. Note that , i.e., a is the distance along L from at which the conditional probability of firing equals one half. For all values of , the estimate of a is close to the distance along L between and the separatrix, which equals 0.05. In other words, the probability of firing, if the path starts at the intersection of L with the separatrix, is about 1 / 2. The estimate of b increases with respect to , and the conditional probability approaches a step function as the amplitude of the noise goes to zero. A step function would correspond to the firing being represented by a first passage time of a fixed threshold.

Fig. 6

Conditional probability of spiking when crossing the line for different values of the noise amplitude . The red dots are individual nonparametric estimates and the blue curve are the fitted curves given by (4.5) (color figure online)

Table 1

Estimates of regression parameters for the conditional probability of firing in the original space (a, b) and in the transformed coordinates based on the additive noise

σ_0	0.001	0.002	0.003	0.004	0.005	0.006	0.007	0.008	0.009	0.01	0.015
a	0.050161	0.050268	0.049946	0.049760	0.049816	0.050001	0.049862	0.049411	0.049078	0.048559	0.046142
b	0.001028	0.002099	0.003192	0.004310	0.005281	0.006459	0.007478	0.008844	0.009877	0.011068	0.017722
a^∗	0.630282	0.631624	0.627576	0.625240	0.625935	0.628262	0.626516	0.620859	0.616673	0.610148	0.579777
b^∗	0.012918	0.026372	0.040106	0.054158	0.066352	0.081158	0.093960	0.111127	0.124107	0.139075	0.222676

To simplify calculations we will work on the transformed coordinates . Then the distance l between (0, l) and (0, 0) in transforms to the distanceand the conditional probability of firing Eq. (4.5) transforms towhere and .

ISI distributions

The comparison of the original stochastic FHN model (3.1) and the two LIF models (4.3) and (4.4) can be performed by studying the ISI statistics. Namely, one first simulates the trajectories of the system (3.1) with starting points close to the fixed point until the first spiking time, and thereafter resets to the starting points. Due to Theorem 3.2, we can simplify the simulation by choosing the starting point at exactly . This was done 1000 times, and the time of the first firing was recorded. A histogram for this data is shown in Fig. 7. The ISI-distribution of Eq. (4.3) is computed as follows (the ISI-distribution of Eq. (4.4) is computed similarly). Let be the first firing time. We computed the density of the distribution of in terms of the conditional hazard rate (Ditlevsen and Greenwood 2013),This function is the density of the conditional probability, given the position on L is r at time t, of a spike occurring in the next small time interval, given that it has not yet occurred.

Fig. 7

The estimated ISI distributions from our approximate LIF models (4.3) and (4.4) with the firing mechanism compare well with the estimated ISI histogram of FHN (3.1) reset to 0 after firings. (color figure online)

Notice that the estimated conditional probability of firing (4.6) is calculated in one cycle of the process, which on average takes time units. Therefore, we estimate the hazard rate asOn the other hand, from standard results from survival analysis, see e.g. Aalen and Borgan (2008) we know that the density of the firing time can be calculated asDue to the law of large numbers, for fixed t, we can numerically determine the density (4.8) up to any desired precision by choosing n and M large enough through the expressionHere ( are M realizations of , and the integral has been approximated by the trapezoidal rule. The results are illustrated in Fig. 7 for , using . The estimated ISI distributions from our approximate LIF models (4.3) and (4.4) with the firing mechanism compare well with the estimated ISI histogram of FHN (3.1) reset to 0 after firings.