Edge detection of remote sensing image based on Grünwald-Letnikov fractional difference and Otsu threshold

1.   Introduction

Among stochastic neuronal models the Leaky Integrate-and-Fire (LIF) model with colored noise (^[1,2]) is a refined version of the classical one driven by white noise (^[3,4]). The standard LIF model cannot explain the high variability in the neuronal response to stimulations and the adaptation phenomena (^[5,6,7,8]), hence more accurate models are needed. The introduction of memory can be achieved in various and different ways. For instance, in ^[9] two different methods of inclusion of memory are presented: Langevin equations with colored noise and fractional Langevin equations. The latter have been also studied in ^[10,11], but we will focus in this work on the first one, which can be used to introduce correlated inputs in the model.

Indeed, correlated inputs have been widely considered to explain adaptation phenomena. One way to introduce correlated inputs is to consider dependence between the excitatory and inhibitory stimuli, as done in ^[12] and ^[13], in which IF and Hodgkin-Huxley models with mutually correlated inputs are considered. Another way to introduce correlated inputs in neuronal models is to consider time-correlated stochastic inputs. This kind of approach is first used in ^[14] and then in ^[1]. Time-correlated inputs are also described in ^[15] with the aid of a binary correlated stochastic process. Finally, let us also recall that time-correlated inputs have been used also to model stochastic currents in spatial models, as done in ^[16].

We mainly refer to the model presented in ^[1], in which time-correlated stochastic inputs are modelled by using a colored noise, that is to say a noise induced by an Ornstein-Uhlenbeck process, in contrast with the white one, which is induced by a Brownian motion. In particular, such colored noise could arise from a diffusive approximation of the synaptic current (see ^[17]). A specific investigation of a colored noise LIF model has been conducted, for instance, in ^[2]. In particular, the processes involved in the model are strictly linked to integrated Gaussian processes. Indeed, we will see that the colored noise induces an integrated Gaussian term in the membrane potential process. This term is non-Markov, then the membrane potential process is not a Gauss-Markov process. However, approximating Gauss-Markov processes can be used to estimate the firing rate of the neuron (see ^[2]). Moreover, the integrated Gaussian term could be further investigated by using some recent techniques on integrated Gaussian processes (see ^[18,19]), which could lead to another estimation of the firing rate of this model. Such investigation will be carried on in future works.

Here, in Section 2, we consider a stochastic model based on a couple of stochastic equations including a correlated input (the colored noise) and a current term. It is derived by applying a linear approximation as in ^[20] to the stochastic model obtained combining those of ^[21] and ^[22]. We give specific justifications about the choice of the involved parameters and, differently from the previous works, we study the model varying the correlation time of the colored noise.

Furthermore, in our investigation we distinguish two different kinds of colored noises, depending on two different reset mechanisms. The first one is linked to the cases in which the noise is induced by some inner mechanism of the neurons (for instance the activation of a gate), and then has to be reset together with the membrane potential process. For this reason we will call it endogenous noise. Moreover, in such case we will consider a non-stationary noise and for modeling it we adopt a non-stationary Ornstein-Uhlenbeck process.

The second one is instead linked with an external noise: it could be seen as the correlated noise of synaptic or injected inputs. For this reason, even if the membrane potential is reset, the inputs are not affected by this reset and so does the noise. Thus we will call it exogenous noise. Finally, in such case, since it does not depend on the dynamics of the neuron, we could assume that the noise is in its stationary state. Then, we will consider a stationary Ornstein-Uhlenbeck (OU) process to model this noise.

We remark that these two cases allow to specialize the proposed model in the attempt to make it more close to some neurophysiological evidences. Moreover, from the mathematical point of view, our aim is also to investigate the model when a non-stationary or a stationary OU process is used to model the noise.

Furthermore, in the model we also include a current input term. We give motivations about the choice of the kind of the current term; then, we show the equivalence between two different ways of including the current term in the model. Indeed, it is possible to consider the current term as an additive term directly in the differential dynamics of the membrane potential process (as in the model here investigated), or, differently, it is possible to consider the current term included in the dynamics of the noise in such a way that an encoded current is in the noise itself (see, for instance, ^[2] and ^[23]). The latter case can be considered to model the dynamics of a neuron embedded in a neuronal network (^[24,25]). Specifically, in Appendix A the link between the model considered in ^[2] and ^[23] and the model here proposed is explained.

In Section 3, for the general case, we give the expressions for the mean and covariance functions of the membrane potential process and for the colored noise process. Then, in Sections 4 and 5 we provide the specialized formulae of the mean and covariance of the membrane potential process in the endogenous and in the exogenous noise case, respectively. By numerical evaluations we plot the surfaces of the covariance functions in order to compare also graphically the results.

Finally, we consider the difficulty that arises in such models related to simulation. Indeed, models involving colored noise also involve differential equations with non-regular terms, hence the classical Euler scheme could not converge. In Section 6, we propose a simulation algorithm to solve this difficulty based on a random Euler scheme given in ^[26] and use it to obtain some numerical approximations of the probability density functions of the first firing time. In particular, we give some measurements of the approximation error in terms of maximum relative error with respect to the mean and the variance function of the process. We also use these measures to compare our algorithm with the classical Euler one. Finally, with these simulations, we will also show in a qualitative way the role of the correlation of the noise in the firing activity.

2.   The main model

We first consider the model based on the following pair of stochastic equations:

where $\eta(t)dt$ is the colored noise: in (2.1) it is in place of the white noise $dW(t)$ as usual in the stochastic differential equation (SDE) of a LIF model. Specifically, the colored noise stands for a correlated input process, $C_m$ is the membrane capacitance, $g_L$ is the leak conductance, $V_L$ is the resting level potential, $I(t)$ is the input current, $W(t)$ is a standard Brownian motion, $\eta_\infty$ is the resting value of $\eta(t)$, $\tau$ is the correlation time of $\eta(t)$ and $\sigma$ is the intensity of the noise. A spike is emitted when $V(t)$ attains a threshold $V_{th}$. The reset mechanism after a spike is fundamental to realize a neuronal model for the firing activity. However, different mechanism of reset of the noise term can be used to describe different physiological phenomena. We will consider:

$(a)$ the case of an endogenous noise, i.e. if $V(t) > V_{th}$ then $V(t^+)\to V_0$ and $\eta(t^+)\to \eta_0$

$(b)$ the case of an exogenous noise, i.e. if $V(t) > V_{th}$ then $V(t^+)\to V_0$ and $\eta(t^+)\to \eta(t) \ \hbox{(no reset in the noise)}.$

Case $(a)$: some modeling justifications.

To better understand the nature of the reset mechanism given in $(a)$, let us consider a conductance based model that takes in consideration the leak, synaptic and muscarinic $K^+$ current. Following the lines of ^[21], in particular Equation $1$ and Section $2.2.2$ for the expression of the muscarinic $K^+$ current, we have

where $p(t)$ is a gating variable and $I_{syn}(t)$ is the synaptic current. For the gating variable $p(t)$, we refer to [22,Equation $11$], thus we have

To obtain a more tractable model, we proceed to linearize the gating variable equation, following the reduction method presented in [20,Section $2.4$ and Appendix $A$]. In particular we have from [20,Equation $37$]

where $\tau_p$ is the (mean) time constant of $p$ and $p_\infty$ is the (mean) equilibrium state of the gating variable, and then, by approximating $V(t)-V_K \simeq \Delta V$ with $\Delta V$ constant, we also have

Moreover, since we are using a fixed threshold model instead of an adaptive threshold one, the gating variable has to be reset together with the membrane potential. In particular if $V(t) > V_{th}$, then $V(t^+) \to V_{reset} = V_0+\delta V$ and $p(t^+) = p_{reset} = p_\infty+\delta p$ where $\delta V$ and $\delta p$ are small perturbations on the initial values that are due to the low spiking activity of a neuron which is subject to a slow muscarinic $K^+$ current (see ^[20]). By simple algebraic transformations, one can recognize in this model Equations (2.1) and (2.2). Indeed, one can set

and, finally we obtain $\eta(t) = g_M p(t)\Delta V$.

Specifications of the parameters as functions of the correlation time of the noise.

In order to specify the above parameters $\eta_\infty, \eta_0, \tau$ in the model (2.1)-(2.2), we need to know $\Delta V, p_\infty$ and $\tau_p$. These quantities can be estimated by considering an average $\overline{V}$ of the process $V$ and using the formulas (see [20,Table 1])

Typically, one could assume $\overline{V} \simeq -60$mV (mean membrane voltage in vivo) and then obtain $\tau_p \simeq 200$ms and $p_\infty \simeq 0.08$. However, we will vary these values, studying the model for different values of $\tau_p$ and determining $\overline{V}$ by inverting formula (2.8) and then $p_\infty$ from formula (2.9).

Moreover, depending on eventual modeling necessity, one could decide to reset also the input current term $I(t)$. Finally, if we assume last spike is far away in time from the initial time of the simulation, we can pose $\delta p \simeq 0$, since $p$ stabilize itself near $p_\infty$ before a spike.

Case $(b)$: some modeling justifications.

Case $(b)$ is quite different. Indeed, it can be used to describe phenomena in which the synaptic current is itself noisy and correlated. Time-correlated synaptic currents have been already considered in LIF models in ^[14] and in ^[1], while they can be used in more complicated models, such spatial models with depolarization (see ^[16]). In particular, colored noise is preferred to white noise in the synaptic currents when the synapses exhibit a finite relaxation time (see, for instance, ^[17], in particular Equations $3.5$ and $3.6$).

Another important difference between the two cases $(a)$ and $(b)$ will be the choice of $\eta_\infty$ and $\eta_0$. In the endogenous noise case $(a)$, the values of $\eta_0$ and $\eta_\infty$ have their effects on the dynamics of the model. Following ^[20], one could choose $\eta_\infty \in \mathbb{R}$, while $\eta_0 = \eta_\infty+\delta \eta$, where $\delta \eta \in \mathbb{R}$ is a perturbation on the resting value $\eta_\infty$. In the exogenous case $(b)$ we do not have any reset mechanism on $\eta$, hence one could consider $\eta_0$ as a free parameter. However, if we want $\eta$ to be a baseline noise, it could be convenient for it to be stationary. To obtain the stationarity of $\eta$, one could consider $\eta_\infty = 0$ and then

(cf. for instance the initial value of Eq.(1.3) in ^[27]).

The current term.

For the input current term $I(t)$, let us consider

This kind of current term mimics the behavior of a synaptic current with exponential decay (see [17,Section $2.2.2$], and [24,25]). Exponential currents have been also used to describe synaptic activity in non-linear models such as the Quadratic Integrate-and-Fire model (see [28]) or the Exponential Integrate-and-Fire model (see [29]). However, our current exhibits a zero baseline justified by the assumption of a zero mean synaptic input (see Equation $3.1$ of ^[17]). Alternatively, such kind of current can be also justified by considering a single synaptic input and an exponential synaptic kernel (see Equation $7.4$ of ^[30]). Moreover, this current is deterministic. Indeed, in the endogenous noise case $(a)$ the noisy part of the current is neglected in favor of the endogenous noise of the gating variable, while in the exogenous noise case $(b)$ the noisy part of the current is explicitly written as the coupled process $\eta(t)$, after using the diffusion approximation given in [17,Equation $3.1$]. One can also use periodic currents instead of exponential ones (see, for instance, [31,32]). However, we will focus on such currents in future works, together with some deeper analysis of the firing times probability densities and the construction of an adaptive threshold model (to substitute the reset mechanism) as done in ^[33].

Here, at first we consider the current in (2.10) as a synaptic current with $I_0 = 0.5$nA, that increases the voltage for a few mV, (see left side of Figure 1 and Figure 2), and $\beta = 5$ms (^[17], ^[30]). Then, we also consider the current in (2.10) as an injected current (as done, for instance, in ^[10]) that can have a greater initial value than the one of a single synaptic input. In particular, to perform simulations of firing times, we set its initial value $I_0 = 3$nA, since a strong excitatory input is required to observe a spike in an admissible time. Indeed, consider that in (2.1) the stochastic input $\eta(t)$ can act like an inhibitory current whose modulus can reach values up to $2.5$nA in its confidence interval of the $99\%$. Furthermore, in (2.10) we set $\beta = 200$ms in order to have a decay time that is comparable to the value of the correlation time $\tau$ ($\tau\sim 200$ms, according to ^[20]).

About the probability density function of firing times, a Gauss-Markov (GM) process can be used to approximate the solution process $V(t)$ (see [2]); hence a Volterra integral equation (see [34]) can be exploited to obtain some approximations of the firing densities.

Here, we determine the mean and covariance functions of the $V(t)$ processes (obtained by solving Equations (2.1) and (2.2)) in order to investigate the different behaviors of the neuronal potential when it is subject to the two different stimulations of cases $(a)$ and $(b)$. Moreover, by simulation procedures we also provide estimates of firing densities in order to give further useful comparisons.

3.   General results

3.1.   The process $V(t)$

Before working with a specific case, let us show what can we say about a generic solution of Eq. (2.1), where $\eta$ is a generic GM process with almost surely continuous paths. We have the following proposition.

Proposition 1. Let $V$ be a strong solution of (2.1) with $V_0 \in \mathbb{R}$, $\eta$ a GM process with almost surely continuous paths adapted to the filtration generated by $\eta_0$ and the Brownian motion $W$, $I(t)$ a continuous function and $\theta: = \frac{C_m}{g_L}$. Hence we have

In particular, we also have

and

Proof. To obtain (3.1), one has just to observe that for fixed $\omega \in \Omega$, Eq. (2.1) is a linear non-autonomous and non-homogeneous ordinary differential equation (ODE), hence it is easy to find a solution (which exists almost surely since $\eta(t)$ has continuous paths for almost every $\omega \in \Omega$, hence $s \mapsto e^{\frac{s}{\theta}}\eta(s)$ is integrable in $[0, t]$).

For Eq. (3.2), we just need to observe that $t \mapsto V_0\, e^{-\frac{t}{\theta}}+V_L\left[1- e^{-\frac{t}{\theta}}\right]$ is a deterministic function, hence it coincides with its own expectation, while

where the second equality follows from Fubini's theorem, since $\eta(s)$ is in $L^2(\mathbb{P})$ for any $s \in [0, t]$, thus $\mathbb{E}[|\eta(s)|] < +\infty$ and, being $\eta$ a GM process, $s\mapsto \mathbb{E}[\eta(s)]$ is continuous, hence integrable in $[0, t]$.

Finally, Eq. (3.3) follows from the fact that

hence

where the last equality follows from Fubini's theorem, since $\eta(u)\in L^2(\mathbb{P})$ for any $u \in [0, \max\{t, s\}]$, hence $\eta(u)\eta(v)\in L^1(\mathbb{P})$ for any $(u, v)\in [0, t]\times [0, s]$, and, being $\eta$ a GM process, $Cov(\eta(u), \eta(v))$ is continuous and then integrable in $[0, t]\times [0, s]$.

From this proposition we also have that the solution process $V$ involves the integral of a GM process. Such integrals are widely studied in ^[18,19].

3.2.   The colored noise process $\eta$

Now let us consider the process $\eta$ solution of (2.2). We know how to write the process $\eta$ in $[0, +\infty)$ as

and then

Hence we know that $\lim_{t \to +\infty}m_\eta(t) = \eta_\infty$ independently from the choice of $\eta_0$. In the endogenous noise case $(a)$, we have that $\eta_0 = \eta_\infty+\delta \eta$, so

and

In the exogenous noise case $(b)$, we pose $\eta_\infty = 0$ and $\eta_0 = \frac{\sigma}{\tau}\int_{-\infty}^0e^{\frac{s}{\tau}}dW(s)$, which is a zero mean Gaussian random variable with variance given by $\frac{\sigma^2}{2\tau}$, that guarantees the stationarity. In particular in this case we have

and $\mathbb{E}[\eta(t)] = 0$. From now on let us denote $m_\eta(t) = \mathbb{E}[\eta(t)]$ and $Cov(\eta(t), \eta(s)) = c_{\eta}(t, s)$. Moreover, since $\eta$ is a GM process, we can define the covariance factors $h_1$ and $h_2$ to be, for $s < t$,

and the ratio that is given by

In the exogenous case, it will be useful the define the function $\rho(t) = r(t)-r(0)$ since $r(0)\not = 0$ in such case.

4.   The model with endogenous noise

Let us consider the endogenous noise case $(a)$. Writing Eq. (2.1) in its integral form, we have

In the endogenous noise case, the covariance function of $\eta$ is given by

Hence we can choose the covariance factors $h_1$ and $h_2$ to be

to obtain a ratio $r$ given by

By using Doob's formula (see ^[35])

for some Brownian motion $\widetilde{W}$ (despite being equal in law, it could be different from the process $W$), we can exploit Eq. (4.1) as

Moreover, recall that in this case the process $\eta$ has to be reset after a spike.

Note that if in Eq. (4.1) $g_L = 0$ and $I(t)\equiv 0, \forall t > 0$, the process $V(t)$ is the integrated over time of a non-stationary Ornstein-Uhlenbeck process. Specifically, in ^[19] a detailed analysis of a similar process (in terms of integrals of GM processes) was done. We can say that those results can be used for the above model with $g_L\equiv 0$ and when the function $I(t)\equiv 0, \forall t > 0$.

Here, by using Proposition 1, we obtain the mean and covariance functions for the special model considered in Eqs. (2.1)-(2.2) with the current $I(t)$ as in Eq. (2.10).

Corollary 1. The stochastic process $V(t)$ of model based on Eqs. (2.1)-(2.2) with $V_0$ real number, $\eta_0 = 0$ and $I(t)$ as in Eq. (2.10) has the following mean function (for $\beta, \tau \not = \theta$)

and covariance function for $t \ge s$

where we set $\theta = \frac{C_m}{g_L}$, $\alpha_1 = \beta-\theta$, $\alpha_2 = \tau-\theta$ and $\alpha_3 = \tau+\theta$.

Proof. These two expressions follow respectively from (3.2) and (3.3) with $m_{\eta}$ and $c_{\eta}$ given in (3.5) and (4.2) and $I$ given in (2.10), after some straightforward but cumbersome calculations.

From the expression of the covariance in (4.4), we can exploit the variance function

where $\theta, \alpha_2, \alpha_3$ are defined in Corollary 1.

Plots of the mean and the variance function are shown in Figure 1. Here, we can see that the mean function admits the horizontal asymptote $\lim_{t \to +\infty}\mathbb{E}[V(t)] = V_L-\eta_\infty\theta/C_m$ (see in the Figure 1, left). Note that the initial excitatory effect of the current $I(t)$ is tempered by the inhibitory one of the endogenous noise $\eta(t)$, being the latter definitely visible for large times. For the variance function, we have $\lim_{t \to +\infty}\mathbb{D}[V(t)] = \frac{\sigma^2\theta^2}{2C_m^2(\tau+\theta)}$, hence its horizontal asymptote decreases as $\tau$ increases, as we can see in the Figure 1 (right). Indeed, as $\tau$ increases, the process has a lower variance, but the covariance decreases more slowly, as we can see in Figure 3.

Also, simulations can be performed. Approximations of the First Passage Time (FPT) probability density function (pdf) were obtained in Section 6. Indeed, the FPT pdf is useful to model the firing density of the neuron.

5.   The model with exogenous noise

In the exogenous noise case $(b)$, the covariance function of $\eta$ is given by

As we specified before, in this case the process $\eta$ is stationary and, as we can see from Eq. (5.1), the function $c_{\eta}(s, t)$ depends only on $|t-s|$. Moreover, the covariance factors can be chosen as

while the ratio $r$ and the function $\rho(t): = r(t)-r(0)$ are given by

Hence, by using the representation of $\eta(t)$ given in Eq. (4.3) and the fact that $\widetilde{W}(r(t))\overset{d}{ = }\widetilde{W}(\rho(t))+\widetilde{W}(r(0))$, we obtain from (4.1):

where $\widetilde{W}(r(0))$ is not almost surely $0$ and it is correlated with $\int _0 ^t h_2 (s) \widetilde{W}(\rho(s)) ds$. Let us also recall that this time, even if we use a different Brownian motion $\widetilde{W}$, we only obtain equality in law.

Here, by still using Proposition 1, we can also obtain the mean and the covariance function of $V$.

Corollary 2. The stochastic process $V(t)$ of model based on Eqs. (2.1)-(2.2) with $V_0$ a real number, $\eta_0 = \frac{\sigma}{\tau}\int_{-\infty}^0e^{\frac{s}{\tau}}dW(s)$, $\eta_\infty = 0$ and $I(t)$ as in Equation (2.10) has the following mean function (for $\beta, \tau \not = \theta$)

and covariance function for $t\ge s$

where we set $\theta = \frac{C_m}{g_L}$, $\alpha_1 = \beta-\theta$, $\alpha_2 = \tau-\theta$ and $\alpha_3 = \tau+\theta$.

Proof. To obtain Eq. (5.2) we just have to use (3.2), recalling that $I$ is given in Eq. (2.10) and $\mathbb{E}[\eta(t)] = 0$ for any $t \ge 0$.

For Eq. (5.3), one just have to use (3.3) where $c_\eta$ is given in (5.1).

From the expression of the covariance in (5.3), we can exploit the variance function

where $\theta, \alpha_2, \alpha_3$ are defined in Corollary 2.

It is important to observe that since $\mathbb{E}[\eta(t)] = 0$, the mean function of $V$ does not depend on $\tau$. Moreover, it admits an horizontal asymptote given by $\lim_{t \to +\infty}\mathbb{E}[V(t)] = V_L$. A plot of the mean function is given in Figure 2 (left). Moreover, since also in this case $\lim_{t \to +\infty}\mathbb{D}[V(t)] = \frac{\sigma^2\theta^2}{2C_m^2(\tau+\theta)}$, the horizontal asymptote of the variance function decreases as $\tau$ increases, as it can be seen in Figure 2 (right). From Figure 4, we can conclude that even in this case, while the variance is lower, the covariance decreases slowly as $\tau$ increases.

Comparing plots in Figure 1 and Figure 2 it is possible to see how the values of $\tau$ affect the behaviour of the variance functions in both cases of noise. Furthermore, in Figure 3 and in Figure 4, as $\tau$ varies, the effects on the surfaces of the covariance functions are also put in evidence.

6.   Simulations and firing densities estimates

Suppose we want to simulate the process $V(t)$ which is part of the solution of the equations

where $I(t)$ is a sufficiently regular function (in our case $I \in C^\infty$). $\eta(t)$ is an Ornstein-Uhlenbeck process, hence its paths are $\beta$-Hölder continuous almost surely for any $\beta < \frac{1}{2}$.

Specifically, the noise $\eta(t)$ is defined on some filtered probability space $(\Omega, \mathcal{F}, \mathcal{F}_t, \mathbb{P})$ hence we have that for $\omega \in \Omega$, $V(\cdot, \omega)$ is the realization of the process $V$ related to the random event $\omega$. If we fix $\omega \in \Omega$ (hence if we consider a single realization of $V$), the equation

is a linear non-homogeneous and non-autonomous Ordinary Differential Equation (ODE) that could be written in the form

Let us then fix

to write the ODE as

However, $f(t, x, \omega)$ is only $\beta$-Hölder continuous in $t$ (except for a null-set with respect to $\omega$) for $\beta < \frac{1}{2}$, hence the classical Euler scheme could not converge. Indeed, to assure convergence of a classical numerical scheme, one needs some regularity assumptions on $f$, not only with respect to $x$ (which are needed to assure the existence of the solution), but also with respect to $t$. In this case, the $\beta$-Hölder continuity of $f$ could be not enough to guarantee the convergence of a classical Euler scheme and then we decide to adopt a Monte-Carlo based scheme. Indeed, since $f$ is a Caratheodory function, we can still use a random Euler scheme (as done in ^[26]). However, to use this scheme we should know exactly the values of $\eta(t, \omega)$. Indeed the random Euler scheme is based on a Monte-Carlo quadrature algorithm (see ^[26] for details), in which the integrand function has to be evaluated in uniformly distributed random times inside the integration domain. Thus, since we need to numerically simulate also $\eta$ (with an Euler scheme as exploited in ^[36]), we have to consider approximatively uniformly distributed random times. The idea of the algorithm is the following: we first simulate $\eta(t)$ with a classical Euler scheme in a finer lattice; then we simulate $V(t)$ by using the random Euler scheme given in ^[26] where the random times for the Monte-Carlo quadrature procedure are chosen on the simulation lattice of $\eta(t)$.

The simulation scheme.

Hence, to simulate the process, let us fix three natural numbers $\delta, \varepsilon, m$ with $\delta < \varepsilon$ and let us define $\Delta t = 10^{-\delta}$ and $\Delta s = 10^{-\varepsilon}$ (in particular $\Delta s$ is a finer time step with respect to $\Delta t$ and $\Delta t/\Delta s \in \mathbb{N}$, while the choice of $m$ is linked to the precision of the Monte-Carlo quadrature scheme). Moreover, let us denote $V_n = V(n\Delta t)$ and $\eta_l = \eta(l\Delta s)$.

Step 1 Initialize $V_0$ and $\eta_0$;

Step 2 Suppose we have simulated $V_n$ and $\eta_l$ such that $n \Delta t = l\Delta s$. We want to simulate $V_{n+1}$;

Step 3 Simulate $\eta_h$ until $h\Delta s = (n+1)\Delta t$: to do that suppose we have already simulated $\eta_{h-1}$ and simulate $\xi \sim \mathcal{N}(0, 1)$, then pose

Step 4 Simulate a vector $U = (U_1, \dots, U_m)$ of i.i.d. uniform random variables in $[0, \Delta t]$;

Step 5 Consider the values $\widetilde{U}_i$ obtained from $U_i$ by rounding up the $\varepsilon$-th decimal digit;

Step 6 Simulate the value $V_{n+1}$ as

where $l(h)$ is defined in such a way that $l(h)\Delta s = n \Delta t+\widetilde{U}_h$ for $h = 1, \dots m$;

Step 7 Return to Step 2.

In particular observe that the quantity $l(h) \in \mathbb{N}$ is well-defined. Indeed, since $\widetilde{U}_h$ is rounded up to the $\varepsilon$-th decimal digit, it can be written as $\widetilde{U}_h = u_h10^{-\varepsilon}$ with $u_h \in \mathbb{N}$. Thus $\Delta t = 10^{-\delta} = 10^{-\varepsilon}10^{\varepsilon-\delta}$ where $\varepsilon-\delta > 0$, and $\Delta s = 10^{-\varepsilon}$; finally $l(h) = n 10^{\varepsilon-\delta}+u_h$.

Finally, let us remark that the reset mechanism can be easily implemented by controlling whenever the simulated process $V_n$ exceeds the threshold $V_{th}$.

Some simulation results

As example, we produced some sample paths of $V$ by using such algorithm in the endogenous noise case, without considering the reset mechanism. These paths are shown in Figure 5, together with a representation of their mean function and a confidence interval determined by using the variance function.

To give an estimate of the error, let us denote with $\overline{V}_n$ the sample mean of a certain number of simulated processes at the time $t_n$ and $D_n$ their sample variance. We will give a measure of the error in terms of maximum relative error, that is to say the quantities.

In particular we will compare our simulation algorithm with the classical Euler scheme, which is given by

where $\xi_n$ is a standard normal random variable. It is interesting to see that since the first step in $V_n$ of Euler method is given by

where $V_0$ and $\eta_0$ are fixed real numbers, $V_1$ is always deterministic, thus $D_1 = 0$ and we have $E_D \ge 1$. As we can see from Table 2, the maximum relative error with respect to the mean $E_M$ seems to be stable at almost $0.1\%$ for our simulation algorithm, while it is almost $1\%$ with the classical Euler scheme. Moreover, it is interesting to see that even for $m = 1$ the maximum relative error with respect to the variance $E_D$ is lower that the one of the Euler scheme. In particular, such error decreases as $m$ increases up to $m\simeq 50$, reaching $E_D \simeq 1.6\%$.

Such simulation algorithm can be useful to estimate the FPT of $V$. Indeed, in Figure 6 we have shown some simulated pdf of these FPT for three different values of $\tau$. In this figure we can see how different values of $\tau$ produce different pdfs. In particular as $\tau$ increases, we can see from this figure that the mode of the FPT shifts to the right. However, a more careful and appropriate analysis of the FPT densities has to be carried on to study how $\tau$ affects such functions.

7.   Conclusions

We have considered a modification of the LIF model with a colored noise which could represent different neuronal mechanisms (such as gating variables or correlated noise for an injected current). This model is described is Section 2 and is constituted of two coupled equations: a linear SDE and a linear ODE with an irregular term. Moreover, in this section we described the different reset mechanisms for the modelling process. In Section 3 we study some general properties of solutions of these coupled equations. In particular, we distinguish between an endogenous noise (in Section 4) and an exogenous noise (in Section 5), and we determine the mean, the covariance and the variance functions of the process $V(t)$ in such cases obtained by adopting a non-stationary and stationary Ornstein-Uhlenbeck process, respectively. Finally, in Section 6 we provide a simulation algorithm for such process, in which we had to keep in mind the irregularity of the terms of the ODE. Such simulation algorithm can be used to obtain some numerical approximation of the FPT of $V(t)$ to a threshold. This last simulation also gives us a validation of the usability of the model. Indeed, the parameters described in Table 1 and $\beta$ are intrinsic of the problem we are studying, while $\sigma$ and $\tau$ should be estimable parameters ($\sigma$ is linked to the variability of the process, hence it can be estimated by using the spread of the pdf of the FPT). The simulations of the FPT cast light on the effect of $\tau$: as $\tau$ grows the mode of the pdf of the FPT shifts to right. Hence $\tau$ is an observable parameter and it could be estimated from the observation of the FPTs.

Acknowledgments

This work is partially supported by GNCS. We acknowledge the constructive criticism of anonymous reviewers on an earlier version of this paper.

Conflict of interest

The authors declare there is no conflict of interest.

A. An equivalent model with the current term included in the colored noise

Another interesting choice for the input current is given by

for some stimuli $s \mapsto S(s)$. In this case, supposing $\eta_\infty = \eta_0 = 0$ we could consider the process $\xi(t) = -\eta(t)+I(t)$ which solves the equation

This kind of stimuli has been considered in ^[23] to describe the case in which the input current received by the neuron has already been codified by other neurons in the network. In particular we have the following equivalence.

Proposition 2. If $I(t)$ is a $C^1$ function then the model given by Equations (2.1)-(2.2) is equivalent to the one given by

that is to say:

$i$ For any $I \in C^1$ there exists a continuous function $t \mapsto S(t)$ such that if $V(t)$ is solution of (2.1) then $V(t)$ is also solution of (A.2);

$ii$ For any $S \in C^0$ there exists a function $I \in C^1$ such that if $V(t)$ is solution of (A.2) then $V(t)$ is also solution of (2.1).

Proof. We gave a closed form for the solution of (2.1) in Equation (3.1). Moreover, we have a closed form for the solution of (A.2) given by

Hence if we want $V$ to be solution of (2.1) and (A.2) simulaneously, we need

that is implied by $\xi(t) = I(t)-\eta(t)$. Moreover, we gave a closed form for a solution of (2.2) in Equation (3.4). We can also give a closed form for a solution of (A.3):

where the $-$ sign can be substituted to the $+$ by symmetry of the Itô integral term. Hence, to obtain the relation $\xi(t) = I(t)-\eta(t)$, one has to pose $\xi_0 = -\eta_0$, and

We have then showed $ii$, since Equation (A.7) gives us the function $I$ we have to choose in order to assure that a solution of (A.2) is also solution of (2.1).

To show $i$, suppose $I \in C^1$ and observe that, from (A.7),

Applying the time derivative operator to this equation we have

thus we finally have

concluding the proof of $i$.