A theoretical connection between the noisy leaky integrate-and-fire and escape rate models: the non-autonomous case

TL;DR

This paper establishes a theoretical link between the noisy leaky integrate-and-fire model and the escape rate model in neuroscience, especially under time-dependent inputs, providing a unified framework that explains their statistical similarities.

Contribution

It introduces a new general stochastic framework that encompasses both models, revealing their connection in non-autonomous conditions.

Findings

The framework unifies the two models under a common theoretical structure.

It explains the observed statistical similarities between the models.

The results have implications for understanding neuronal noise processes.

Abstract

One of the most important challenges in mathematical neuroscience is to properly illustrate the stochastic nature of neurons. Among different approaches, the noisy leaky integrate-and-fire and the escape rate models are probably the most popular. These two models are usually chosen to express different noise action over the neural cell. In this paper we investigate the link between the two formalisms in the case of a neuron subject to a time dependent input. To this aim, we introduce a new general stochastic framework. As we shall prove, our general framework entails the two already existing ones. Our result has theoretical implications since it offers a general view upon the two stochastic processes mostly used in neuroscience, upon the way they can be linked, and explain their observed statistical similarity.

Equations

$$\frac{\text{mbox}\mathrm{d}}{\text{mbox}\mathrm{d}t}v(t)=(\mu (t)-v)+\sigma \xi (t),$$

$$⟨\xi (t)⟩=0,⟨\xi (t)\xi (t^{\prime })⟩=\delta (t-t^{\prime }).$$

$$\text{mbox}ifv\ge 1\text{mbox}thenv\mapsto v_r.$$

$$\frac{\partial }{\partial t}p(t,v)+\frac{\partial }{\partial v}[(\mu (t)-v)p(t,v)]-\frac{{\sigma }^2}{2}\frac{{\partial }^2}{\partial v^2}p(t,v)=\delta (v-v_r)r(t),$$

$$p(t,1)=0,$$

$$r(t)=-\frac{{\sigma }^2}{2}\frac{\partial }{\partial v}p(t,1).$$

$$\underset{v\to -\infty }{\lim }[(v-\mu (t))p(t,v)+\frac{{\sigma }^2}{2}\frac{\partial }{\partial v}p(t,v)]=0.$$

$$p(0,v)=p_0(v),$$

$${\int }_{-\infty }^1{p}_0(v)\text{mbox}\mathrm{d}v=1.$$

$${\int }_{-\infty }^{1}p(t,v)\text{mbox}\mathrm{d}v=1.$$

$$\frac{\text{mbox}\mathrm{d}}{\text{mbox}\mathrm{d}t}a(t)=1.$$

$$\mbox{ If the spike occurs at time $t$ then }a\mapsto 0.$$

$$\frac{\partial }{\partial t}n(t,a)+\frac{\partial }{\partial a}n(t,a)+S(t,a)n(t,a)=0,$$

$$n(t,0)=r(t)={\int }_0^{\infty }S(t,a)n(t,a)\text{mbox}\mathrm{d}a,$$

$$n(0,a)=n_0(a).$$

$${\int }_0^{\infty }n(t,a)\text{mbox}\mathrm{d}a=1,$$

$${\int }_0^{\infty }{n}_0(a)\text{mbox}\mathrm{d}a=1.$$

$$\left\{\begin{array}[]{ll}\frac{\mbox{\rm\,d}}{\mbox{\rm\,d}t}v(t)=(\mu(t)-v(t))+\sigma\xi(t),\\ \frac{\mbox{\rm\,d}}{\mbox{\rm\,d}t}a(t)=1,\end{array}\right.$$

$$\text{mbox}ifv\ge 1\text{mbox}thenv\mapsto v_r\text{mbox}anda\mapsto 0.$$

$$⟨\xi (t)⟩=0,⟨\xi (t)\xi (t^{\prime })⟩=\delta (t-t^{\prime }).$$

$$\frac{\partial }{\partial t}\pi (t,a,v)+\stackrel{\text{Drift}}{\overbrace{\frac{\partial }{\partial a}\pi (t,a,v)+\frac{\partial }{\partial v}[(\mu (t)-v)\pi (t,a,v)]}}-\stackrel{\text{Diffusion}}{\overbrace{\frac{{\sigma }^2}{2}\frac{{\partial }^2}{\partial v^2}\pi (t,a,v)}}=0.$$

$$\pi (t,a,1)=0,$$

$$\underset{v\to -\infty }{\lim }(-\mu +v)\pi (t,a,v)+\frac{{\sigma }^2}{2}\frac{\partial }{\partial v}\pi (t,a,v)=0.$$

$$\underset{a\to +\infty }{\lim }\pi (t,a,v)=0.$$

$$\rho (t,a)=-\frac{{\sigma }^2}{2}\frac{\partial }{\partial v}\pi (t,a,1),$$

$$\pi (t,0,v)=\delta (v-v_r){\int }_0^{\infty }\rho (t,a)\text{mbox}\mathrm{d}a:=\delta (v-v_r)r(t).$$

$$\pi (0,a,v)={\pi }_0(a,v).$$

$${\int }_{-\infty }^1{\int }_0^{\infty }{\pi }_0(a,v)\text{mbox}\mathrm{d}a\text{mbox}\mathrm{d}v=1,$$

$${\int }_{-\infty }^1{\int }_0^{\infty }\pi (t,a,v)\text{mbox}\mathrm{d}a\text{mbox}\mathrm{d}v=1.$$

$$p(t,v)={\int }_0^{\infty }\pi (t,a,v)\text{mbox}\mathrm{d}a,$$

Full text

Keywords: Neural noise; Noisy Leaky Integrate-and-Fire model; Escape rate

1 Introduction

The noisy aspect of neuronal behavior [15] makes the use of stochastic processes unavoidable in mathematical neuroscience. A generally accepted formalism to express the variability present in neural dynamics has not been found yet, but among the most notable models dealing with this aspect stand out the noisy leaky integrate-and-fire (NLIF) [7, 8], and the escape rate models [17]. Of course, the use of each of these models comes with its own advantages.

The NLIF - a model going back to Lapique’s work in 1907 [3, 1, 6] - describes neurons as simple electrical circuits consisting in a capacitor in parallel with a resistor driven by a noisy input. Written in the language of stochastic differential equations, it gives rise to a Langevin equation [20] completed with a reset mechanism to mimic the onset of an action potential. Related to a stochastic variable satisfying the Langevin equation, stands the probability density function (pdf) that expresses the likelihood that the state variable takes on a specific value [16]. For the NLIF model, its associated pdf follows the well-known Fokker-Planck (FP) equation [16].

The escape rate models - a formalism proposed by W. Gerstner and J.L. van Hemmen in 1992 [17] - introduce a rather different approach, in which the neuron’s dynamics remain deterministic but the generation of an action potential occurs randomly. The firing probability is then described by a stochastic intensity of firing or hazard function which depends on the momentary distance between the membrane potential and the formal firing threshold [17, 18]. In this setting, the neuron can fire even without reaching the formal spiking threshold or remain quiescent after crossing it [18]. The corresponding pdf to the escape rate model satisfies a pure transport equation along the deterministic trajectories, and a term to express the probabilistic nature of spike initiation is added. The so-obtained pdf follows an age structured (AS) system, where the age variable accounts for the time elapsed since the neuron had its last spike [17].

At a first sight, the two stochastic processes reminded above seem essentially different, both in form and in modeling assumptions. Despite their contrasting nature, it has been observed that they can behave in a surprisingly comparable fashion [24]. The main goal of our paper is therefore to investigate the existence of a common ground between the two formalisms.

In our previous work [13, 14], we made a first step in this direction. We have highlighted an unforseen relationship between the two modeling frameworks, but, we have left several open questions, especially in the context of time dependent stimuli. To push further our previous studies, we will now introduce a new stochastic model that describes the evolution of the membrane’s potential of a neuron and its time passed by since the last firing event. We will show that the system satisfied by the associated pdf generates the two previously introduced formalisms: FP equation and AS system.

The paper is structured as follows: We shall first remind the two above mentioned models as they are usually encountered in the literature. Next, we introduce a two dimensional stochastic description of the neural state, and, starting from there, we show that our general stochastic framework entails the two already existing ones. The corresponding stationary system is discussed next and a special solution of the non stationary system is proposed later on. We finish with a discussion on the biological and theoretical implications of our findings.

2 Two standard models for neural noise

2.1 Fokker-Planck equation

The first process -the NLIF- that we shall refer to is modeled as a Langevin equation which describes the membrane potential of a neuron in the subthreshold regime. The subthreshold dynamics are described by:

[TABLE]

where $\mu(t)$ is the mean of the external stimuli, $\sigma$ is the noise intensity of the external stimuli, and $\xi$ is a Gaussian white noise

[TABLE]

The firing mechanism of the neuron is expressed by a discontinuous reset process which says that when the membrane potential reaches the threshold value, the potential is reset to a fix given value:

[TABLE]

By associating to the stochastic variable $v$ a pdf $p(t,v)$, it is found that the pdf satisfies the FP equation [brunel2015fokker]:

[TABLE]

where an absorbing boundary condition at the spiking threshold is imposed in order to express the immediate reset of the variable,

[TABLE]

and the flux at this boundary gives the firing rate:

[TABLE]

At the lower boundary, the no-flux condition is imposed:

[TABLE]

An initial repartition is assumed as known:

[TABLE]

for which the normalization condition is imposed, due to the interpretation of $p$ as a pdf:

[TABLE]

It can be easily shown that, once the condition (7) is assumed, then the conservation property of the solution to the system (3)-(6) takes place at any time:

[TABLE]

In Fig. 1, a simulation of the NLIF model (1), (2) is presented. The first panel (Fig. 1A) corresponds to the simulation of the same neuron under different noise realizations. As expected, the corresponding dynamics are different and the neurons spike at different times. In the second panel (Fig. 1B), the spiking activity of the same cell is extracted from a bigger number of trials. Finally, in the last panel (Fig. 1C) the corresponding spiking activity is compared with the firing rate obtained from the FP model (3)-(6).

In Fig. 2, a comparison is made under a stimulus that is time dependent. The first panel (Fig. 2A) shows the stimulus’ time evolution. In the second panel (Fig. 2B) the spiking activity of the cell is displayed, from which we can extract the firing rate. A comparison with the FP equation is shown in the last panel (Fig. 2C).

2.2 Age-structure formalism

The escape rate models consider a deterministic trajectory combined with a probabilistic firing mechanism [17]. Implicitly, the state variable will depend only on the last firing time. Without restraining the generality, we can consider the state variable for the deterministic trajectories as the time elapsed since the last spike, that we will generically call age:

[TABLE]

The probability that an action potential occurs during a small time interval is computed via the hazard function $S(t,a)$ that gives the probability of spike of a neuron that has at time t the age $a$. More precisely, during the time interval $(t,t+dt)$, a spike occurs with probability $S(t,a(t))dt$ and, consequently, the age is reset to zero immediately after:

[TABLE]

The associated probability density function, $n(t,a)$, satisfies the following PDE [17, 18]:

[TABLE]

where the first two terms give the pure transport along the deterministic trajectories (9) and the last term accounts for the probability of spiking.

The reset mechanism is modeled as a non-local boundary condition:

[TABLE]

where the right hand side of the above equation defines the firing rate of the system.

As before, an initial repartition is taken as known:

[TABLE]

The conservation property of the system (11)-(13) takes place as well. Namely,

[TABLE]

for any $t$, as soon as

[TABLE]

We present in Fig. 3 a simulation of the model (11)-(13). Again, we have made a comparison between the stochastic process (red curve) given by (9), (10) and the evolution in time of the density function (black curve) given by the system (11)-(13). The first panel (Fig. 3A) corresponds to the simulation of the same neuron. As expected, due to the probabilistic firing process, the corresponding dynamics have different spike times. In the second panel (Fig. 3B), the spiking activity of the same cell is extracted from a large number of trials. In the last panel (Fig. 3C) the corresponding spiking activity is compared with the firing rate given by the AS system (12).

In Fig. 4, a numerical simulation of the stochastic process defined by (9), (10) for the same age-dependent death rate is illustrated. Note that the neuron never fires exactly at the same age, since its probability to fire (escape) is purely stochastic. The first panel (Fig. 4A) shows the time evolution of the stimulus. In the second panel (Fig. 4B), the spiking activity of the cell is shown, from which the firing rate is extracted. A comparison between the firing rate directly computed and the one extracted from the AS system is presented in the last panel (Fig. 4C).

3 An age and potential structured model. Connections with FP and AS systems

3.1 The age-and-potential structured model

We shall introduce below a more general process which can be seen as a generalization of the two reminded formalisms. Let us consider a stochastic process that models both variables, i.e. spike times and membrane potential. More precisely, let us consider that, between two consecutive firing times, the potential variable evolves according to the subthreshold dynamics given by the NLIF model (1) and the age of the neuron grows linearly with time according to the dynamics given by (9). We suppose that the neuron will fire whenever a given potential threshold value is reached. Whenever this happens, both state variables are reset: the potential is reset to $v_{r}$ and the age to zero. We deal therefore with the following two-dimensional stochastic process:

[TABLE]

completed by the firing mechanism:

[TABLE]

As in the previous section, $\mu(t)$ is the mean of the external stimulus, $\sigma$ is the noise intensity of the external stimulus, and $\xi$ is a again Gaussian white noise

[TABLE]

The corresponding probability density function satisfies a two-dimensional Fokker-Planck equation (see [16]):

[TABLE]

For the right choice of the boundary conditions, we turn again to the stochastic system (14), (15). Namely, corresponding to the firing mechanism, we should impose an absorbing boundary condition at the firing threshold:

[TABLE]

and a reflecting boundary at the boundary $v=-\infty$:

[TABLE]

We will also assume that

[TABLE]

The flux at the firing threshold is given by

[TABLE]

and, since once firing occurs both membrane potential and the age of a neuron are reset in the reset value $(0,v_{r})$, the following singular source term arises:

[TABLE]

Finally, an initial distribution is assumed as known:

[TABLE]

The model is now complete, and one can check directly by integration over the state space that, if the initial distribution satisfies

[TABLE]

then, for any $t>0$:

[TABLE]

In Fig. 5, a simulation of the age-potential pdf is shown. The different panels (A-B-C) show the evolution in time of the density. The simulation starts with a Gaussian as initial condition (the first panel of Fig. 5). Under the influence of the drift term, the density function advances in age, which is clearly seen in the plots of Fig. 5. After the spiking process, the age of the neuron is reset to zero and its membrane potential is reset to $v_{r}$. This effect is well perceived in the second panel of Fig. 5.

3.2 Theoretical connections

We are ready to prove now the main result of our paper which shows that both models reminded in the previous sections, i.e. FP model (3)-(6) respectively the AS model (11)-(13), can be obtained from our two-dimensional system. We stress that the AS system derived from the age-and-potential system (16)-(19) has a special age-dependent death rate, as it can be seen below.

3.2.1 Derivation of the Fokker-Planck model

The following result takes place:

Theorem 1

Let $\pi$ a solution to the system (16)-(19). Then, the probability density function defined as

[TABLE]

is the solution to the FP system (3)-(6) with the initial repartition given by

[TABLE]

Remark. Here and below, the solutions are defined in distributional sense. We did choose though to keep the computations at a formal level so that to render the reading easy.

• Proof

The proof is simple and straightforward. Let us integrate the equation (16) with respect to $a$ over the age interval:

[TABLE]

which implies, by using the boundary condition at $a\longrightarrow\infty$:

[TABLE]

Now, using the boundary condition (18) and denoting by

[TABLE]

it follows that, indeed, $p$ is solution to

[TABLE]

where

[TABLE]

The boundary condition is obtained by:

[TABLE]

and the initial condition is defined as

[TABLE]

We end the proof by noticing that, since the conservation property of the probability density $\pi(t,a,v)$ is assumed, we get that

[TABLE]

as soon as

[TABLE]

is fulfilled. The proof is now complete.

3.2.2 Derivation of the age-structured model

The following result takes place:

Theorem 2

Let $\pi$ a solution to the system (16)-(19). Then, the probability density function defined as

[TABLE]

is the solution to the AS system (11)-(13) with the initial repartition given by

[TABLE]

• Proof

By integrating (16) with respect to $v$ over the respective domain, one gets:

[TABLE]

where, as defined above,

[TABLE]

The conservation property imposed on the system in $\pi$ will necessary lead, when integrating (22) with respect to $a$ over the age domain, to:

[TABLE]

As before, the initial condition is defined as

[TABLE]

The system (22)-(24) is an equivalent form of the system introduced in the first section. To see this, there should be reminded the significance of the terms in the initial system.

We start by reminding that $S(t,a)$ was defined as the decay rate of the survivor function, which, in turn was defined as the integral over whole potential space of a neuron that had the last spike at $t-a$; notice that, in terms of our system, this is nothing else than $n(t,a)$. We therefore can force the notation into our system, and write down:

[TABLE]

to arrive to the popularized version of the equation that gives the evolution with respect to the time elapsed since the last spike:

[TABLE]

where

[TABLE]

The proof is now complete.

3.3 The stationary system

Considering the probability that a neuron reaches a potential value $v$ at time $a$, given that at the initial time was in the state $v_{r}$, this probability is shown to satisfy the following problem [25, 26]:

[TABLE]

with the initial condition given by

[TABLE]

The model is completed by an absorbing boundary condition at the firing threshold

[TABLE]

and a reflecting boundary condition at the lower boundary of the domain:

[TABLE]

This problem, known as the first passage time problem, can be therefore thought to express the neuronal probability density function’s dynamics on an inter-spike interval. It is due to this reason that we have kept here the notation $a$ for the time variable, since, in the context introduced above, this time variable is interpreted as the time elapsed since the last spike. All these considerations are made under the assumption that the stimulus $\mu$ is time-independent. As it is well known, the flux at the firing threshold

[TABLE]

has the interpretation of the inter-spike distribution density function, and the corresponding survivor function is given by:

[TABLE]

Denoting now by $\bar{\pi}$ the solution to the corresponding stationary system to the model (16)-(19), i.e.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

we obtain that it is given by

[TABLE]

with

[TABLE]

the stationary firing rate.

It is easy to see that, by integrating (37) with respect to $v$ and $a$, respectively, on the corresponding state spaces, we obtain that the distributions defined by:

[TABLE]

respectively,

[TABLE]

satisfy the corresponding stationary systems to AS and FP models.

In the case of the FP and AS models (3)-(6), respectively, (11)-(13), with a time-independent stimulus $\mu$, we have found in [13, 14] a way to map the solutions one-to-another via integral transforms.

In this particular case, the following result takes place:

Theorem 3

Let us consider the system (16)-(19) with a constant in time stimulus $\mu(t)\equiv\mu$. Then the system admits a separable solution

[TABLE]

where $\varphi$ is the solution to (27)-(30) and $n(t,a)$ is the solution to

[TABLE]

with

[TABLE]

Here,

[TABLE]

with $ISI(a)$ and $P(a)$ defined by (31), (32) respectively.

It is easy to check that the function defined by (38) is indeed a solution to the system (16)-(19). Moreover, by integrating this function with respect to $a$ over the respective state space, we recover the transform of the solution to AS model into the solution of the FP model in case of time independent stimuli [13]. The recovery of the solution to the AS system in this case is trivial, reducing it by integration of (38) with respect to v over the respective state space, to the the tautology $n=n$, by taking into account that the survivor function $P$ is given by (32).

3.4 A special solution for the non-autonomous case

In fact, the result from the previous subsection can be adapted for the time dependent stimulus case, given that a corresponding first passage time problem is correctly defined.

Let us define, therefore, the probability density for a neuron to be at time $t$ in the state $(a,v)$, given that at the firing time $t-a$ it was in the state $(0,v_{r})$ and denote it by $\varphi(t,a,v)$. Then, the corresponding FP equation is given again by (16), but with a different boundary condition:

[TABLE]

The rest of the boundary conditions are the same as above:

[TABLE]

[TABLE]

As initial condition, we consider

[TABLE]

with $\varphi_{0}(a,v)$ solution to the first passage time problem for the constant in time stimulus ($\mu\equiv\mu(0)$) (27)-(30).

The survivor function is then by definition:

[TABLE]

and the corresponding ISI function is given by

[TABLE]

It checks out immediately that

[TABLE]

which is natural given the interpretation of the survivor function. A corresponding relation between the survivor function and the ISI function takes place, namely,

[TABLE]

where $\frac{\rm{D}}{\rm{Dt}}$ is the total derivative with respect to $t$. The hazard rate is then defined as

[TABLE]

Let us consider now the system (25)-26), this time with the hazard rate defined by (45). Then, the following result holds:

Theorem 4

The function defined as

[TABLE]

is a solution to the system (16)-(18), where $n(t,a)$ is a solution to (25)-(26).

• Proof

The proof is, again, straightforward. Let us formally check that equation (16) is satisfied:

[TABLE]

The boundary condition at $a=0$ checks out immediately:

[TABLE]

where

[TABLE]

Note that Theorem 4 does allow an actual integral transform of the solution to the AS system into the solution of the FP equation, similar to those obtained in [13, 14]:

[TABLE]

if we assume that initial repartition is given by

[TABLE]

with $\varphi_{0}$ defined above, so that:

[TABLE]

4 Discussion

The present paper is thought as a continuation of our previous work [13, 14], where we have uncovered a connection between the NLIF and the escape rate models in the context of a time independent stimulus. We proved there that the corresponding solutions to the FP equation and the AS system can be mapped one-to-another via integral transforms. In the present paper, we wish to go further by considering stimuli that are time dependent. Toward this end, we introduced a new stochastic model and we showed that starting from it, it is possible to recover both the NLIF and the escape rate models.

The importance of our finding is mostly theoretical. Our conclusion is that it is always possible to adjust parameters such that the NLIF and the escape rate models to exhibit the same statistics. This finding enables us to find one possible explanation for the reported similarities between the two frameworks [24]. The choice of using any of it would then be influenced only by which of the variables - age or membrane potential - one wants to take into account.

Although the FP equation in neuronal dynamics context has been studied over the past few years and qualitative results regarding its solutions have been proven [9, 10], the AS systems theory has received a lot of attention in the last decades, see monographs [19, 27] for example, and different associated control problems have been considered [2]. There is, therefore, a big advantage from a mathematical point of view to privilege its use in neuronal modeling.

Let us finally conclude by saying that our work opens several pathways for future research. From our side, the most exciting one is the investigation of emergent properties of neural networks. Both formalisms (FP and AS) have been employed to describe neural circuits in the mean-field approximation (see [4, 5, 9] for the FP and see [21, 11, 22, 23] for the AS), and both of them have generated key insights in neuroscience, especially in the quest of the underlying mechanism of neural synchronization [12, 23]. It would therefore be relevant to use our new framework to investigate the links and differences between the two approaches, especially when conclusions are contradictory.