Evolution of moments and correlations in nonrenewal escape-time processes

The theoretical description of nonrenewal stochastic systems is a challenge. Analytical results are often not available or can be obtained only under strong conditions, limiting their applicability. Also, numerical results have mostly been obtained by ad hoc Monte Carlo simulations, which are usually computationally expensive when a high degree of accuracy is needed. To gain quantitative insight into these systems under general conditions, we here introduce a numerical iterated first-passage time approach based on solving the time-dependent Fokker-Planck equation (FPE) to describe the statistics of nonrenewal stochastic systems. We illustrate the approach using spike-triggered neuronal adaptation in the leaky and perfect integrate-and-fire model, respectively. The transition to stationarity of first-passage time moments and their sequential correlations occur on a nontrivial time scale that depends on all system parameters. Surprisingly this is so for both single exponential and scale-free power-law adaptation. The method works beyond the small noise and time-scale separation approximations. It shows excellent agreement with direct Monte Carlo simulations, which allow for the computation of transient and stationary distributions. We compare different methods to compute the evolution of the moments and serial correlation coefficients (SCCs) and discuss the challenge of reliably computing the SCCs, which we find to be very sensitive to numerical inaccuracies for both the leaky and perfect integrate-and-fire models. In conclusion, our methods provide a general picture of nonrenewal dynamics in a wide range of stochastic systems exhibiting short- and long-range correlations.

Figure 1

Sample paths of the model for a power-law adaptation current given by Eq. (3). Top: $X\left(t\right)$ [Eq. (1)], the horizontal dashed lines are at $X=0$ and $x_{\mathrm{th}}=1$. Bottom: $s\left(t\right)$ [Eq. (2)] When $X$ reaches $x_{\mathrm{th}},s$ undergoes a jump of size $\kappa$. The subsequent ISIs $T_k$ have distributions ${\mathcal{F}}_k\left(t\right)$, and the starting values $s_0^{\left(k\right)}$ have distributions ${\mathcal{G}}_k\left(s\right)$ for $k\ge 1$. Parameter values are $\alpha =3.0,\gamma =1.0,\sigma =0.8,I_0=4.0,\kappa =3.0$.

Figure 2

Top: Evolution of the rate $r_k=\frac{1}{{\tau }_k^1}$ (left) and standard deviation given by Eq. (6) (right) of ISIs as a function of $k$. Empty circles: Plain MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Filled circles: moments obtained from numerical solution of FPE using a CN time-stepping scheme. $M={10}^6$ independent MC realizations for each value of $k$. Power-law adaptation [Eq. (3)] with $\alpha =5.5,I_0=6.0,\sigma =1.3,\gamma =1.0,\kappa =5.5$. Bottom: Relative disagreements defined by Eq. (13), where an MC-GS algorithm was used.

Figure 3

Top: Evolution of the rate $r_k=\frac{1}{{\tau }_k^1}$ (left) and standard deviation given by Eq. (6) (right) of ISIs as a function of $k$. Empty circles: Plain MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Filled circles: moments obtained from numerical solution of FPE using an Euler time-stepping scheme. $M={10}^6$ independent MC realizations for each value of $k$. Single exponential adaptation [Eq. (4)] with ${\tau }_a=1.0,I_0=5.0,\sigma =1.0,\gamma =1.0,\kappa =1.0$. Bottom: Relative disagreements defined by Eq. (13), where an MC-GS algorithm was used.

Figure 4

PDF ${\mathcal{F}}_k$ of the $k\mathrm{th}$ ISI $T_k$ [Eq. (5)]. The symbols are MC simulations ($M={10}^6$ MC realizations) as indicated in the legends. Solid black lines: PDF obtained by numerical solution of the FPE, Eq. (8). In the left panel, the distributions are practically indistinguishable after $k=2$. Top: Power-law adaptation [Eq. (3)]. Bottom: Single exponential adaptation [Eq. (4)]. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation. Both panels show results for plain MC simulations.

Figure 5

PDF ${\mathcal{G}}_k$ of $s_0^{\left(k\right)}$ [Eq. (7)]. The symbols are MC simulations ($M={10}^6$ MC realizations) as indicated in the legends. The black solid line is ${\mathcal{G}}_1$ obtained from Eq. (11) or ${\mathcal{G}}_2$ (${\mathcal{G}}_3$) obtained using Eq. (18). In the left panel, the distributions are practically indistinguishable after $k=2$. Top: Power-law adaptation [Eq. (3)]. Bottom: Single exponential adaptation [Eq. (4)]. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation. Both panels show results for plain MC simulations.

Figure 6

Conditional FPT densities $\mathcal{H}(\lambda ,\nu )$ given by Eq. (17) relating the initial value of the adaptation current to the distribution of the following ISI. Computed using numerical solutions to the FPE Eq. (8). Top: Power-law adaptation, computed using a CN time-stepping scheme. Bottom: Exponential adaptation, computed using an Euler time-stepping scheme. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation.

Figure 7

MC simulations with GS boundary correction for the joint density of $s_0^{\left(n\right)}$ (i.e., $n\mathrm{th}$ peak value of adaptation current following the $n\mathrm{th}$ FPT) and $T_n$ for different values of $n$. Top: Power-law adaptation, $n=2$. Bottom: Exponential adaptation, $n=3$. $M={10}^6$ MC realizations. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation.

Figure 8

Top: Expectation $\mathbb{E}\left(T_n{T}_{n+1}\right)$ of adjacent ISIs for different values of $n$. Left: power-law adaptation. Right: exponential adaptation. $M={10}^6$ MC realizations. Empty circles: Plain MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Pentagons: Eq. (24). For power-law adaptation, the pentagons are on top of the empty triangles and filled circles and hence not visible. Filled circles: Eq. (27). Diamonds: Eq. (28). The diamonds are nearly on top of the filled circles and hence not visible. The PDE results were obtained using a CN scheme (power-law adaptation) or an Euler time-stepping scheme (exponential adaptation). The vertical error bars show the MC error for a $>99.99\%$ confidence interval. Bottom: Relative disagreements defined by Eq. (13), where for power-law adaptation, the GS boundary corrected MC algorithm was used, and for exponential adaptation, the plain MC algorithm was used. For the iFPT quantity, Eq. (24) was used. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation.

Figure 9

Top: ${\mathcal{Q}}_1(n,1)$ [defined in Eq. (22)] of adjacent ISIs for different values of $n$. Left: power-law adaptation. Right: exponential adaptation. Empty circles: MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Filled circles: Eq. (22), CN time-stepping scheme for power-law adaptation and Euler time-stepping scheme for exponential adaptation. Bottom: Relative disagreement defined by Eq. (13), a plain MC algorithm was used to obtain the relative disagreement. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation.

Figure 10

Top: ${\mathcal{Q}}_2(n,1)$ [defined in Eq. (23)] for adjacent ISIs for different values of $n$. Left: power-law adaptation. Right: exponential adaptation. Empty circles: MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Filled circles: Eq. (23), CN time-stepping scheme for power-law adaptation and Euler time-stepping scheme for exponential adaptation. Bottom: Relative disagreement defined by Eq. (13), a plain MC algorithm was used to obtain the relative disagreement. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation. Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation.

Figure 11

Serial correlation coefficient at lag $n=1$ [defined in Eq. (21)] for adjacent ISIs for different values of $n$ for the LIF model [Eq. (1)]. Left: power-law adaptation. Right: exponential adaptation. Empty circles: Plain MC simulations of Eqs. (1) and (2). Triangles: MC simulations with GS boundary correction. Filled circles: PDE results, SCC computed using Eq. (27). Pentagons: PDE results, SCC computed using Eq. (24). CN time-stepping scheme for power-law adaptation and Euler time-stepping scheme for exponential adaptation. Diamonds: PDE results, SCC computed using Eq. (28). Parameter values as in Fig. 2 for power-law adaptation and as in Fig. 3 for exponential adaptation. Even if the joint expectations shown in Fig. 8 agree well, this does not imply that the SCC will be well approximated; the small correlation values [i.e., the difference between $\mathbb{E}\left(T_n{T}_{n+1}\right)$ and ${\mathcal{Q}}_1(n,1)$] for the two examples lead to large discrepancies in the SCCs, which are especially significant for the case of power-law adaptation.

Figure 12

Serial correlation coefficient at lag $n=1$ [defined in Eq. (21)] for different values of $n$ for the PIF model [Eq. (30)]. The dashed horizontal line is the stationary SCC given by Eq. (31). Empty circles: Plain MC simulations of Eq. (30). Pentagons: Eq. (24). Filled circles: Eq. (27). The PDE results were obtained using a CN time-stepping scheme. $M={10}^6$ MC realizations. Time step $h={10}^{-3}$. Parameter values: $D=0.1,{\tau }_a=5.0,\tilde{\mathrm{\Delta }}=10,I_0=5.5,s_0^{\left(0\right)}=5.0.$