Damped oscillations of the probability of random events followed by absolute refractory period: exact analytical results

TL;DR

This paper provides exact analytical formulas describing damped oscillations in the probability of random events following an absolute refractory period, using a stochastic neuron model as an example.

Contribution

It introduces a novel, simplified analytical approach to describe transient oscillations in stochastic processes with refractory periods, avoiding renewal theory.

Findings

Derived explicit formulas for oscillation amplitude damping.

Demonstrated the approach with a stochastic neuron model.

Provided insights into transient dynamics of refractory processes.

Abstract

There are numerous examples of natural and artificial processes that represent stochastic sequences of events followed by an absolute refractory period during which the occurrence of a subsequent event is impossible. In the simplest case of a generalized Bernoulli scheme for uniform random events followed by the absolute refractory period, the event probability as a function of time can exhibit damped transient oscillations. Using stochastically-spiking point neuron as a model example, we present an exact and compact analytical description for the oscillations without invoking the standard renewal theory. The resulting formulas stand out for their relative simplicity, allowing one to analytically obtain the amplitude damping of the 2nd and 3rd peaks of the event probability.

Equations

$$\nu ={\bar{p}}_s/△t=p_s/(△t+p_s{\tau }_{ref}).$$

$$P_k=p_s\left[1-P(A_{k-n_{ref}}+A_{k-n_{ref}+1}+...+A_{k-1})\right].$$

$$P(A_{k-n_{ref}}+A_{k-n_{ref}+1}+...+A_{k-1})=\sum _{j=k-n_{ref}}^{k-1}P(A_j),$$

$$P_k=p_s(1-\sum _{j=k-n_{ref}}^{k-1}{P}_j)=p_s\left[1-(S_{k-1}-S_{k-n_{ref}-1})\right],$$

$$S_m=\{\begin{array}{l}0,m\le 0,\\ \sum _{j=1}^m{P}_j,m>0.\end{array}$$

$$S_k=P_k+S_{k-1}=p_s(1+S_{k-n_{ref}-1})+(1-p_s)S_{k-1}.$$

$$S_k=S_1{r}^{k-1}+\frac{r^{k-1}-1}{r-1}d=1-(1-p_s{)}^k.$$

$$P_k\equiv P_k^{(1)}=p_s(1-S_{k-1})=p_s(1-p_s{)}^{k-1}.$$

$$P_{\infty }=p_s/(1+n_{ref}{p}_s).$$

$$P_k=\{\begin{array}{l}{p}_s(1-p_s{)}^{k-1},k\le n_{ref}+1,\\ p_s{P}_{k-n_{ref}-1}+(1-p_s)P_{k-1},k>n_{ref}+1.\end{array}$$

$$P_k^{(2)}=(k-n_{ref}-1)p_s^2(1-p_s{)}^{k-n_{ref}-2}+P_k^{(1)}=P_k^{(1)}[1+(k-n_{ref}-1)q],$$

$$P_k^{(3)}$$

$$=P_k^{(1)}[1+(k-n_{ref}-1)q+\frac{1}{2}(k-2n_{ref}-2)(k-2n_{ref}-1)q^2].$$

$$P_k^{(m+1)}=P_k^{(1)}[\frac{1}{m!}\prod _{j=1}^m(k-m\cdot n_{ref}-j)]q^m+P_k^{(m)}=P_k^{(1)}{C}_{k-m\cdot n_{ref}-1}^m{q}^m+P_k^{(m)}.$$

$$\frac{1}{m!}\prod _{j=1}^m(n-j)=C_{n-1}^m,$$

$$P_k=P_k^{(1)}(1+\sum _{i=1}^m{C}_{k-i\cdot n_{ref}-1}^i{q}^i),$$

$$P_k=P_k^{(1)}+\sum _{j=1}^{k-1}{P}_{j}P(A_k|A_j)=P_k^{(1)}+\sum _{j=1}^{k-n_{ref}-1}{P}_j{P}_{k-j-n_{ref}}^{(1)}.$$

$$P_k=P_k^{(1)}+\sum _{i=1}^{k-n_{ref}-1}{P}_{k-i-n_{ref}}{P}_i^{(1)}.$$

$$k_{\max 2}=n_{ref}+1+R,$$

$$P_{\max }^{(2)}\equiv P_{k_{\max 2}}^{(2)}=P_{k_{\max 2}}^{(1)}\frac{q}{u}=\frac{p_s^2}{u}(1-p_s{)}^{R-1}.$$

$$D_{i+1}=P_{\max }^{(i+1)}/P_{\max }^{(i)}.$$

$$D_2=\frac{p_s}{u}(1-p_s{)}^{R-1}.$$

$$x{k}^2-yk-z=0,$$

$$z=-(2n_{ref}+\frac{3}{2})q^2-(n_{ref}+1)(2n_{ref}+1)q^{2}u+q+(n_{ref}+1)qu-u.$$

$$k_{\max 3}=\frac{y+\sqrt{y^2+4xz}}{2x}=2(n_{ref}+1)+R+X,$$

$$X=-\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{1}{u^2}-\frac{(2n_{ref}+1)}{q}-\frac{1}{q^2}}.$$

$$P_{k_{\max 3}}^{(3)}=P_{k_{\max 3}}^{(1)}[1+(n_{ref}+1+R+X)q+\frac{1}{2}(R+X)(1+R+X)q^2].$$

$$D_3=P_{k_{\max 3}}^{(3)}/P_{k_{\max 2}}^{(2)}=p_s(\frac{1}{2}+X+\frac{1}{u})(1-p_s{)}^X.$$

$$P_k^{(rd)}=h(k△t)△t.$$

$$h(t)=\sum _{n=1}^{\infty }{f}_n(t),$$

$$f_n(t)=\underset{0}{\overset{+\infty }{\int }}{f}_{n-1}(x)f(t-x)dx.$$

Full text

There are numerous examples of natural and artificial processes that represent stochastic sequences of events followed by an absolute refractory period during which the occurrence of a subsequent event is impossible. In the simplest case of a generalized Bernoulli scheme for uniform random events followed by the absolute refractory period, the event probability as a function of time can exhibit damped transient oscillations. Using stochastically-spiking point neuron as a model example, we present an exact and compact analytical description for the oscillations without invoking the standard renewal theory. The resulting formulas stand out for their relative simplicity, allowing one to analytically obtain the amplitude damping of the 2nd and 3rd peaks of the event probability.

Keywords: Renewal point process, Event probability, Absolute refractory period, Damped oscillations, Stochastically spiking neuron

1. Introduction

Natural and technical events are often followed by a refractory period: after the event has happened, repeating similar events for some time is either unlikely (relative refractory period) or impossible (absolute refractory period). A characteristic natural example is the refractory nature of the neuron ability to generate electrical pulses, "spikes". In turn, a typical technical example is the so-called dead time for some type of detectors, particularly photodetectors: the dead time is a fixed period after the detector triggering during which it becomes inoperative. Both the neuronal refractoriness and the detector dead time make an essential impact on the statistical distribution of counted events. Assuming that instant identical events occur randomly and that each event is followed by the same absolute refractory period $\tau_{ref}$, it is intuitively clear that average interval $T$ between the events should be a sum, $T=\tau_{ref}+T_{0}$, where $T_{0}$ is the average interval if the refractoriness would vanish. A more nontrivial consequence of the refractory period is damped oscillations of the event occurrence probability as a function of time Perk67 ; WCBJ72 ; Kin06 . These are strongly pronounced at $T_{0}\lesssim\tau_{ref}$. Such transient oscillations arise due to (i) a specific initial condition and (ii) the fact that, in the presence of refractoriness, the random instant events ordered in time are no longer independent from each other. In particular, the probability of a subsequent event depends on the time elapsed from the preceding one. Random point processes of this kind are called renewal processes and are a subject of study of the renewal theory Cox62 . This theory has a method of getting an explicit analytical description for the damped oscillations of the event probability by finding the so-called renewal density Perk67 ; Cox62 . For instance, it has been done repeatedly for the classic example of the Poisson process modulated by the dead time or absolute refractory period Malm47 ; Ric66 ; Muller73 ; Muller74 ; Cant75 ; John83 ; John86 ; Pom99 ; Picib08 ; PRE2010 ; JCN2012 ; NC2018 . However, the resulting analytical formula for the time-dependent probability of the event is quite bulky, and that complicates substantially getting further analytical findings.

In this paper, for the Bernoulli process modulated by absolute refractoriness, we give a compact analytical description of the damped oscillations without invoking the renewal theory. The description is presented in four equivalent forms: three variants of a recurrence formula and the explicit formula as a series. These are quantitatively consistent with numerical simulations and known results of the renewal theory. One of the recurrence formula variants is particularly simple, enabling accurate analytical calculation of the damping coefficients. Notably, these are quite robust against changing the values of the model parameters.

2. Model description and problem statement

For certainty, we consider a model neuron that stochastically emits spikes, each of which is followed by the absolute refractory period. In particular, the neuron can spontaneously emit a spike with probability $p_{s}$ per unit time (i.e., in a given elementary time interval $\triangle t$) so that after the spike emission the neuron becomes temporarily inactive: it cannot emit spikes during the refractory period $\tau_{ref}=n_{ref}\triangle t$, where $n_{ref}$ is a positive integer. The mean occurrence rate of spiking events in the absence of refractoriness is given by $\nu_{0}=1/T_{0}=p_{s}/\triangle t$. At $\tau_{ref}\neq 0$, from the equality $T=\tau_{ref}+T_{0}$, for the mean rate $\nu=1/T$ one gets $\nu=\nu_{0}/(1+\tau_{ref}\nu_{0})$. In fact, this mean rate is equal to the ratio of the total number of spikes to the observation time $T_{obs}$, given that $T_{obs}>>\tau_{ref}$, $T_{0}$. It is also useful to introduce an asymptotic value for the average probability of spike generation at each step, $\bar{p}_{s}=\nu\triangle t=p_{s}/(1+n_{ref}p_{s})$, such that, by the analogy with the formula for $\nu_{0}$,

[TABLE]

The algorithm for simulating the neuron’s dynamics is quite simple and as follows. Dividing the observation interval $T_{obs}$ by $N$ equal steps $\triangle t$, $T_{obs}=N\triangle t$, these time steps are numbered by a sequence of natural numbers starting with 1. At each step, a random number $\xi$, uniformly distributed from zero to one, is generated and compared with the given probability $p_{s}$ of generating a spike. If $\xi\leq p_{s}$, it is assumed that spike has been generated at this time step. After the event, the neuron cannot emit a next spike during refractory period $\tau_{ref}=n_{ref}\triangle t$. The event of spike generation at the arbitrary $k$-th step is further denoted by $A_{k}$. For definiteness, the initial state of the neuron is chosen as a moment when the neuron has just left the refractory state after emitting a spike.

Performing either a large number of repeated passes of the observation interval $T_{obs}$ for a single neuron or a single pass for the large ensemble of independent neurons, one gets the statistical distribution of the event occurrence in the entire sequence of $N$ time intervals. Normalizing this distribution by the number of either the passes or the neurons in the ensemble, one can obtain the probability distribution $P_{k}$ of spike generation at $k$-th elementary time step of the observation interval. Due to the refractory period, the event probability, as a function of time, exhibits damped oscillations with the average period equal to $\tau_{ref}$ (Fig. 1).

In turn, the analytical problem consists in finding the probability $P_{k}$ of spike generation at each $k$-th elementary time step so that in the asymptotic limit $k\rightarrow\infty$ one would obtain $P_{k}\rightarrow\bar{p}_{s}$.

3. Derivation of the exact analytical formula for $P_{k}$

The probability of spike generation at the $k$-th time step ($k=\overline{1,N}$), $P_{k}\equiv P(A_{k})$, is equal to the product of $p_{s}$ and the probability that in the interval from $k-n_{ref}$ to $k-1$ inclusively no spike has been emitted,

[TABLE]

The events of spike generation in $n_{ref}$ consecutive time intervals are pairwise incompatible. Therefore the probability of the sum in Eq. (2) equals the sum of probabilities

[TABLE]

and the sought-for probability at the $k$-th step is determined in a recurrent manner, with the recursion period equal to the refractory period (cf. JNM01 ),

[TABLE]

where, by definition,

[TABLE]

Such a definition of $S_{m}$ allows us to directly generalize the formula for $P_{k}$ to the range $1\leq k\leq n_{ref}$. The closed formula for $S_{k}$ has the form of a linear recurrent sequence

[TABLE]

The resulting formula (4) accurately describes the numerical statistics (Fig. 2).

One should note three important consequences. First, for the range $1\leq k\leq n_{ref}+1$, where $S_{k}=p_{s}+(1-p_{s})S_{k-1}$, the sum $S_{k}$ can be easily found explicitly, as it is an arithmetic-geometric progression of the form $S_{k}=rS_{k-1}+d$ with $S_{1}=p_{s}$, $r=1-p_{s}$, and $d=p_{s}$. According to the formula for the explicit form of the $k$-th term of this progression,

[TABLE]

Hence, from Eq. (4) the probability of spike generation in the interval $1\leq k\leq n_{ref}+1$ is given by

[TABLE]

Here and below, the upper index $m$ in parentheses in $P_{k}^{(m)}$ indicates the sequential number of the current refractory interval counted from the initial moment $t=0$. The probability (7) refers to the geometric distribution and, given the initial condition, naturally coincides with the probability of the first spike generation at the arbitrary $k$-th time step, such that ${\displaystyle\sum\limits_{k=1}^{\infty}}P_{k}^{(1)}=1$.

The second consequence is that formula (4) makes it easy to obtain the asymptotic probability value $P_{k}$ at $k\rightarrow\infty$. Indeed, the difference $S_{k-1}-S_{k-n_{ref}-1}$ contains $n_{ref}$ terms. At $k\rightarrow\infty$ the probability at the $k$-th step remains practically unchanged. Denoting it as $P_{\infty}$, from Eq. (4) one gets $P_{\infty}=p_{s}(1-n_{ref}P_{\infty})$, whence

[TABLE]

Finally, the third consequence: calculating the adjacent terms $P_{k+1}$ or $P_{k-1}$, alike to Eq. (6), one can exclude sums (5) from Eq. (4) and obtain linear recurrent sequence for $P_{k}$:

[TABLE]

This formula is completely equivalent to Eq. (4). Using Eq. (9), one can easily find the expression for $P_{k}$ in an explicit form within the intervals of $k$ multiples of $n_{ref}+1$.

For example, at $n_{ref}+1\leq k\leq 2(n_{ref}+1)$ one gets

[TABLE]

where $q=p_{s}(1-p_{s})^{-(n_{ref}+1)}$.

Next, at $2(n_{ref}+1)\leq k\leq 3(n_{ref}+1)$ one gets

[TABLE]

It is worth noting that in Eqs. (10) and (12) the numerical coefficient in the highest-order term with respect to $q$ is the so-called triangular number $j(j+1)/2$ at $j=k-2n_{ref}-2$. The inset in Fig. 2 shows the plots for $P_{k}^{(1)}$, $P_{k}^{(2)}$, and $P_{k}^{(3)}$.

By the induction method, one can obtain a formula for $P_{k}^{(m+1)}$, which is valid within the range $m(n_{ref}+1)\leq k\leq(m+1)(n_{ref}+1)$:

[TABLE]

Here, equality

[TABLE]

where $n>m\geq 1$, and $C_{n}^{k}=\frac{n!}{k!(n-k)!}=\dbinom{n}{k}$ is the standard binomial coefficient, has been used to get the last formula in Eq. (13).

In turn, Eq. (13) allows one to derive the general explicit expression for $P_{k}$ in the polynomial form:

[TABLE]

where $m$ is the integer part of the rational number $k/(n_{ref}+1)$.

Another equivalent formula for $P_{k}$ can be obtained as follows. Denote $P(A_{k}|A_{j})$ the conditional probability of spike generation at the $k$-th step, provided that the previous spike was generated at the $j$-th step. At times greater than the refractory period, i.e. at $k>n_{ref}$, the probability of generating a subsequent spike depends only on the moment of the previous spike. Therefore, taking into account the initial condition, $P(A_{k}|A_{j})=P_{k-j-n_{ref}}^{(1)}$ at $k-j>n_{ref}$, and $P(A_{k}|A_{j})=0$ at $k-j\leq n_{ref}$.

If $k>n_{ref}+1$, probability $P_{k}$ can be written as a sum of two terms: (i) the probability $P_{k}^{(1)}$ that a spike will be emitted for the first time on the $k$-th step and (ii) the probability that at least one spike has been emitted previously. The latter has the form of a convolution and follows from the total probability formula. In particular, one gets

[TABLE]

Notably, using substitution $i=k-j-n_{ref}$, one can virtually swap the indices of the multipliers under the sign of the sum in Eq. (16), while the formula does not change its numerical value:

[TABLE]

Despite the different appearance in relation to Eqs. (4) and (9), the last two formulas lead to the same numerical results. In fact, these can be derived from the recurrent formula (9): e.g., Eq. (17) can be obtained by successively substituting values $P_{k-1}$, $P_{k-2}$, $\ldots$ , $P_{k-n_{ref}-1}$ into Eq. (9) and using the definition (7) for $P_{k}^{(1)}$.

4. Damping of the oscillations

The explicit formulas (10) and (12) allow analytical calculation of the relative damping of the second and third peaks of $P_{k}$. Finding the location of these peaks corresponds to solving a linear and quadratic algebraic equation, respectively. In particular, for $(n_{ref}+1)\leq k\leq 2(n_{ref}+1)$, assuming that $k$ is a continuous variable and calculating $dP_{k}^{(2)}/dk=0$ with the use of Eq. (10), one gets the location of the second peak,

[TABLE]

where $R=1/u-1/q$, $u=\ln[1/(1-p_{s})]$, and $q=p_{s}(1-p_{s})^{-(n_{ref}+1)}$. Substituting $k_{\max 2}$ into Eq. (10), we get the amplitude of the second peak

[TABLE]

The damping can be defined as the ratio of amplitudes of the adjacent decreasing peaks,

[TABLE]

Given that $k_{\max 1}=1$ and $P_{k_{\max 1}}^{(1)}=p_{s}$, for the second peak we get

[TABLE]

For the third refractory interval, $2(n_{ref}+1)\leq k\leq 3(n_{ref}+1)$, calculation of $dP_{k}^{(3)}/dk=0$ using Eq. (12) results in quadratic equation,

[TABLE]

where $x=\frac{1}{2}q^{2}u$, $y=q^{2}+(2n_{ref}+\frac{3}{2})q^{2}u-qu$, and

[TABLE]

A suitable solution of Eq. (22) is the root

[TABLE]

where

[TABLE]

Substituting $k_{\max 3}$ into Eq. (12) for $P_{k}^{(3)}$ yields

[TABLE]

After elementary but cumbersome calculations (given in the Supplementary Material), one can get a compact analytical expression for the damping coefficient of the third peak,

[TABLE]

Numerical calculations have confirmed the validity of the obtained formulas. Below we have also listed the numerical values of the relevant quantities, calculated by the above formulas, for two examples shown in Fig. 1 and Fig. 2.

For the first example, at $\tau_{ref}=2$ ms ($n_{ref}=\tau_{ref}/\triangle t=200$ at $\triangle t=0.01$ ms) and $p_{s}=0.1$, we get $q\approx 1.6\times 10^{8}$, $u\approx 0.1$, $R\approx 9.5$, $k_{\max 2}=210$, $P_{k_{\max 2}}^{(2)}\approx 0.04$, $k_{\max 3}=420$, $P_{k_{\max 3}}^{(3)}\approx 0.03$, $D_{2}\approx 0.39$, and $D_{3}\approx 0.74$.

For the second example, at $\tau_{ref}=5$ ms ($n_{ref}=500$) and $p_{s}=0.01$, we get $q\approx 1.5$, $u\approx 0.01$, $R\approx 99$, $k_{\max 2}=599$, $P_{k_{\max 2}}^{(2)}\approx 0.004$, $k_{\max 3}=1196$, $P_{k_{\max 3}}^{(3)}\approx 0.003$, $D_{2}\approx 0.37$, and $D_{3}\approx 0.75$.

It is seen that the numerical values $D_{2}\approx 0.4$ and $D_{3}\approx 0.75$ are fairly robust against changing the parameters. Thus, the magnitude of the second peak is approximately equal to 40% of the magnitude of the first. In turn, the value of the third peak is approximately 75% of the value of the second or 30% of the value of the first peak.

5. Comparison with the renewal theory

In the framework of the renewal theory Cox62 ; Perk67 , when the time intervals between events are independent random variables, one can standardly find the time dependence of the event probability, provided that an event occurred at the initial moment. In particular, the sought-for probability is expressed through the so-called renewal density $h(t)$,

[TABLE]

In turn, the renewal density $h(t)$ is determined by the distribution density $f(\tau)$ for the time interval $\tau$ between the successive events Cox62 ; Perk67 ,

[TABLE]

where $f_{1}(t)\equiv f(t)$ and for $n\geq 2$ functions $f_{n}(t)$ are given by the recursive convolution

[TABLE]

Qualitatively, functions $f_{n}(t)$ are the distribution densities of so-called $n$-th order time intervals between events Perk67 : the first-order interval is the time elapsed from some benchmark event to the first following event, the second-order interval is the time elapsed from the benchmark event to the second following event, etc. The $n$-th order interval is therefore a sum of $n$ successive first-order intervals and is spanned by $(n+1)$ successive events.

At asymptotically large time $t\rightarrow+\infty$, $h(t)$ is saturated Cox62 ,

[TABLE]

Here $\rho$ is the mean rate of events, defined as the inverse mean interval between the successive events,

[TABLE]

There are many, likely independent, examples of applying the renewal theory results to the case, where the event occurrence is the Poisson process modulated by the constant dead time (in our notations, absolute refractory period $\tau_{ref}$) or, equivalently, the distribution density of time intervals between the successive events has a form of the shifted exponential distribution,

[TABLE]

where $\nu_{0}=p_{s}/\triangle t$ and $\theta(\ldots)$ is the unit step function: $\theta(x)=1$ at $x\geq 0$ and $\theta(x)=0$ otherwise. In particular, to the best of our knowledge, the first relevant paper dates back to 1947 Malm47 and has been followed by both in-depth mathematical studies Cox62 ; Ric66 ; Muller73 ; Muller74 ; Cant75 ; Pom99 ; Picib08 ; PRE2010 and applied studies for neuroscience Perk67 ; John83 ; JCN2012 ; NC2018 ; TMC78 ; JASA85 ; John86 ; Koch93 ; Berry98 ; Reich98 ; Kas01 (see also overviews in John96 and Gerst02 ). Below, we briefly outline and compare the previous results with our findings.

For $f(t)$ given by (32), using the Laplace transform, one can reduce Eq. (29) to the following expression

[TABLE]

which is the probability density function for the Erlang (or gamma) distribution, enabling computation of the renewal density $h(t)$ and the sought-for probability (27). The plot of the function $P^{(rd)}_{k}$ given by Eq. (27) with

[TABLE]

where $n_{\max}=10^{3}$ and $f_{n}(t)$ is determined by (33), is shown in Fig. 3. The time dependence $P^{(rd)}_{k}$, accurate to an offset equal to the refractory period, completely coincides with $P_{k}$ calculated by mutually equivalent formulas (4), (9), (15), (16), (17) in Section 3. The offset arises due to the different initial conditions: an exit from the refractory period at $t=0$ in our model and an event occurrence at $t=0$ in the standard renewal-density approach.

Finally, using (30), (31) and (32), one can get the asymptotic value of $h(t)$,

[TABLE]

which is consistent with formulas (1) and (8) for $\nu$ and $P_{\infty}$, respectively.

6. Discussion

In most cases of multi-unit systems, such as biological neuronal networks with neurons as the units, the absolute refractory period $\tau_{ref}$ is not exactly the same for each unit. For instance, it may be randomly distributed so that each neuron has its own value of $\tau_{ref}$ (e.g., see Fig. 2a in Avi13 ). Alternatively, $\tau_{ref}$ may be updated randomly after every event Pet21 . Then increasing the variance of the refractory period distribution would apparently blur the damped oscillations. Nevertheless, one can safely set the absolute refractory period as a positive constant value, same for all neurons, in the majority of neuronal network models without losing their main predictions. Therefore the obtained analytical results may be directly used as a simple test of correctness for the neuronal network models with stochastic spontaneous activity of single neurons, provided that the stochasticity mimics internal neuronal dynamics rather than incoming synaptic currents from some "external" neurons (e.g., see PRL93 ; PRE95 ; PRE98 , cf. Newh10 ). Indeed, in the former case, one can approximate the probability of spontaneous spike generation as nearly time-independent, constant quantity. On the contrary, in the latter case, this probability may be essentially time-dependent, even for cortical neurons in the high-input regime, when they receive hundreds of excitatory synaptic inputs during each interspike interval SN1998 ; NRN2003 ; Kuhn2004 . Nevertheless, the obtained results may be still applicable for the stochastic firing-rate-based neuron models (e.g., JN2008 ; SRM1 ; SRM2 ), where the neuronal dynamics are entirely described by a time-dependent probability of spike generation. It should also be mentioned that there exists a different type of stochasticity induced by the multilayer network structure Perc16 ; Perc17 and the results obtained in this paper are not directly applicable to this type.

The test mentioned above is as follows: in the limiting case where all synaptic interactions between neurons are set to zero, one should increase the value of probability of spontaneous spiking by a single neuron (implying that this value has been made the same for all neurons of the network) until the pseudo-network activity (or, more precisely, the fraction of neurons emitted spikes during the elementary time step of simulation dynamics) would exhibit pronounced damped oscillations. Next, turning to the analytical formulas, one should either compute the time dependence of spiking probability (e.g., using Eq. (9)), and match this with the simulation result or, even simpler, compare the relative values of the first, second, and third peaks of the normalized collective spiking activity of neurons with the values obtained by formulas (21) and (26). If the neuronal network model has been designed and implemented correctly, one gets the complete correspondence, as disconnecting neurons of the network turns it into a statistical ensemble of independent neurons (cf. Truc2017 ). In fact, the plots in Fig. 1 have been originated from two actual realizations of the test applied to the spiking neuronal network model reported in ZP2021 .

In addition, the analytical results may be of interest for developing and testing (i) population-level "neural mass" models and (ii) network-based "cellular automata" models.

The top-down population models, such as the classical Wilson-Cowan model WCBJ72 , usually include the absolute refractory period as a phenomenological parameter. Notably, the damped oscillations due to the absolute refractory period were the first consequence from the original Wilson-Cowan equations compared to their temporally coarse-grained version (see Fig. 3 in WCBJ72 ), though after some reasoning the authors concluded that "the damped oscillation is not likely to have any functional significance" (WCBJ72 , page 9). In a more recent bottom-up approach to neuronal population modeling, sometimes refereed as the refractory density approach Chizh1 ; Chizh2 ; Chizh3 (cf. Omur2000 ; Gerst2000 ; Mat2002 ; Deco2008 ; Non2010 ; Gerst17 ; Dum17 ; PRE20 ), enabling to extract the population dynamics starting with relatively complex, biophysical Hodgkin-Huxley-type neuron models, an occurrence of the damped oscillations due to the refractoriness is the foundation for a standard test on the transient response of the population firing rate to a step-like increase of the common input stimulation.

In turn, a wide spectrum of network-based cellular automata models, which typically constitute randomly-connected simple excitable units with nonlinear interaction, has been elaborated for large-scale numerical and/or analytically-tractable studies of critical collective behavior via avalanches, emergent synchronization and global oscillations of neuronal activity JN2008 ; SRM1 ; SRM2 ; SR2017 (some earlier findings were reviewed in Lind2004 ). In addition to the foregoing models, some relevant models, where interference between the absolute refractory period and stochastic event occurrence may be essential, are the two-state unit model Esc2012 , the three-state unit model PA2003 ; Wood1 ; Wood2 (see also JSM2011 ; Lima19 ), the DeVille-Peskin multi-state model BMB2008 , the cortical branching model Widom2014 , and others Golt10 ; Lee2014 ; NJP15 .

Finally, the analytical results obtained in this paper might be useful for a qualitative analysis of statistical data on both spontaneous and stimulated spiking activity of real biological neurons. For instance, the neurons could have receptive fields subject to an intense step-like stimulus (for the fly’s photoreceptors, such a case was considered in Song17 , see Figs. 5 and 6 there). The references to relevant experimental data are as follows: Figs. 6 and 15 in Pog64 , Fig. 4 in Rod67 , Fig. 1 in Gray67 , Figs. 1, 3, 4, 5, and 6 in Moor70 , Fig. 3 in Rob87 , Fig. 2B in Nin95 , Figs. 2H and 8E in Kita04 , Figs. 4A and 4B in Kita06 , Fig. 2A in Beat12 , Fig. 4A in Lee14 , and Fig. 4A in Amit17 . In turn, direct quantitative comparison of the analytical results with the above experimental findings is hindered by such additional factors as (i) the relative refractory period, which may be time-dependent, and (ii) uncontrollable non-random interaction of the target neurons with surrounding cells Kara2000 ; Kass04 ; Amar06 ; Naw07 ; Maim09 ; EBR11 ; Coh11 . These and other factors Truc05 may have a dominant influence on the heights, widths, and locations of experimental autocorrelogram peaks.

7. Conclusion

The model considered in this paper is the Bernoulli scheme supplemented by the condition of absolute refractoriness followed after each random event. It formally refers to the renewal theory. However, being quite simple, the model allows obtaining useful results without invoking this formalism. In particular, four equivalent analytical descriptions of the damped oscillations of the event probability have been given: (i) recurrent formula (4) through the difference of two sums, (ii) closed recurrent formula (9), (iii) explicit formula (15) in the form of a polynomial, and (iv) recurrent convolution-type formulas (16), (17). It has also been shown that for the Poisson approximation these results accurately coincide with that of the renewal theory. Finally, using the closed recurrent formula, the relative damping coefficients for the second and third peaks of the event probability have been found in the exact analytical form.

It is worth emphasizing that the example of stochastic neuron spiking serves only as a baseline illustration of the obtained mathematical results, which have potentially much wider range of applicability. In general, the damped oscillations studied in this paper are well-pronounced if the absolute refractory period is a dominantly large timescale in the system’s dynamics responsible for the event occurrence.

Data and code availability

The Supplementary Material to this paper contains ready-to-use MATLAB/Octave codes for performing simulations and plotting graphs of the obtained analytical formulas. It also contains a detailed derivation of the damping coefficients for the 2nd and 3rd peaks of the event probability (see Section 4).

Author contributions

A.V. Paraskevov conceived the study, designed the mathematical model and program codes, performed simulations and analytical calculations, prepared graphs, and wrote the manuscript. A.S. Minkin derived the formulas (16) and (17), independently reproduced all the rest results, and verified the manuscript.

Declaration of competing interests

The authors declare no competing interests.

Acknowledgments

A.V. Paraskevov thanks E.Z. Meilikhov, L. Logiaco, and D. Festa for the stimulating discussions.