Stochastic bursting in unidirectionally delay-coupled noisy excitable systems

TL;DR

This paper demonstrates that stochastic bursting occurs in delay-coupled noisy excitable systems due to the interplay of time delays and noise, with analytical descriptions matching simulations.

Contribution

It introduces an analytical framework for understanding stochastic bursting in delay-coupled excitable systems, highlighting the role of timescale separation and leader-follower dynamics.

Findings

Stochastic bursting observed in delay-coupled noisy systems.

Analytical spike statistics and correlations derived.

Good agreement between theory and simulations.

Abstract

We show that \emph{stochastic bursting} is observed in a ring of unidirectional delay-coupled noisy excitable systems, thanks to the combinational action of time-delayed coupling and noise. Under the approximation of timescale separation, i.e., when the time delays in each connection are much larger than the characteristic duration of the spikes, the observed rather coherent spike pattern can be described by idealized coupled point processes with a leader-follower relationship. We derive analytically the statistics of the spikes in each unit, pairwise correlations between any two units, and the spectrum of the total output from the network. Theory is in a good agreement with the simulations with a network of theta-neurons.

Equations

$$\dot{{\theta }_i}=a_i+\cos {\theta }_i+{ϵ}_i(a_{i-1}+\cos {\theta }_{i-1}(t-{\hat{\tau }}_{i-1}))+\sqrt{D_i}{\xi }_i(t).$$

$$H(t)=a+\cos {\Theta }_{sp}(t),$$

$${\Theta }_{sp}(t)=2\arctan \left(\sqrt{\frac{1+a}{1-a}}\tanh \left(\frac{\sqrt{1-a^2}}{2}(t-t_0)\right)\right).$$

$$\frac{\partial P(\theta ,t)}{\partial t}=-\frac{\partial }{\partial \theta }\left[(a+\cos \theta )P(\theta ,t)\right]+D\frac{{\partial }^{2}P(\theta ,t)}{\partial {\theta }^2}.$$

$$P_{st}(\theta )=C{\int }_{\theta }^{\theta +2\pi }\frac{d\psi }{D}{e}^{-{\int }_{\theta }^{\psi }\frac{a+\cos \phi }{D}d\phi }.$$

$$\lambda =C\left(1-e^{-{\int }_0^{2\pi }\frac{a+\cos \theta }{D}d\theta }\right).$$

$$p(ϵ)=\frac{N(ϵ)-N(0)}{N(ϵ)}.$$

$${\mu }_i$$

$$S_i(T)$$

$$Q_i(T)=\{\begin{array}{ll}1-e^{-{\mu }_{i}T},& T<\tilde{\tau },\\ 1-(1-p_i{p}_j)e^{-{\mu }_i\tilde{\tau }-({\lambda }_i+{\lambda }_j{p}_j)(T-\tilde{\tau })},& T\ge \tilde{\tau },\end{array}$$

$$P_{ij}(t,t+\Delta t;t+s,t+s+\Delta t)=\{\begin{array}{ll}{W}_i(t,t+\Delta t)P_{ij}(t+s|t,\Delta t),& s\ge 0,\\ W_j(t,t+\Delta t)P_{ji}(t|t+s,\Delta t),& s<0.\end{array}$$

$$P_{ii}(t+s|t,\Delta t)$$

$$=\sum _{k=0}^{\infty }(p_i{p}_j{)}^k\delta (s-k\tilde{\tau })\Delta t,s\ge 0;$$

$$P_{ii}(t|t+s,\Delta t)$$

$$C_{ii}(s)$$

$$=\frac{1}{T}{\int }_0^{T}dt\underset{\Delta t\to 0}{\lim }\frac{P_{ii}(t,t+\Delta t;t+s,t+s+\Delta t)}{\Delta t^2}={\mu }_i\sum _{n=-\infty }^{\infty }(p_i{p}_j{)}^{|n|}\delta (s-n\tilde{\tau }),$$

$$\begin{array}{r}{S}_{ii}(\omega )=\underset{-\infty }{\overset{\infty }{\int }}{C}_{ii}(s)e^{-i\omega s}ds=\frac{({\lambda }_i+{\lambda }_j{p}_j)(1+p_i{p}_j)}{1+(p_i{p}_j{)}^2-2p_i{p}_j\cos \omega ({\tau }_i+{\tau }_j)}.\end{array}$$

$$P_{ij}(t+s|t,\Delta t)=$$

$$=\sum _{k=0}^{\infty }{p}_i(p_i{p}_j{)}^k\delta (s-k\tilde{\tau }-{\tau }_i)\Delta t,s\ge 0;$$

$$P_{ji}(t|t+s,\Delta t)=$$

$$=\sum _{k=0}^{\infty }{p}_j(p_i{p}_j{)}^k\delta (s+k\tilde{\tau }+{\tau }_j)\Delta t,s<0.$$

$$C_{ij}(s)=$$

$$={\mu }_i\sum _{n=0}^{\infty }{p}_i(p_i{p}_j{)}^n\delta (s-n\tilde{\tau }-{\tau }_i)+{\mu }_j\sum _{n=0}^{\infty }{p}_j(p_i{p}_j{)}^n\delta (s+n\tilde{\tau }+{\tau }_j).$$

$$\begin{gathered} S_{ij}(\omega )=\underset{-\infty }{\overset{\infty }{\int }}{C}_{ij}(s)e^{-i\omega s}ds \\ ={\mu }_i{p}_i\frac{e^{-i\omega {\tau }_i}}{1-p_i{p}_j{e}^{-i\omega \tilde{\tau }}}+{\mu }_j{p}_j\frac{e^{i\omega {\tau }_j}}{1-p_i{p}_j{e}^{i\omega \tilde{\tau }}}. \end{gathered}$$

$$\Re \{S_{ij}\}=\frac{p_i({\mu }_i-{\mu }_j{p}_i{p}_j)\cos \omega {\tau }_i+p_j({\mu }_j-{\mu }_i{p}_i{p}_j)\cos \omega {\tau }_j}{1+(p_i{p}_j{)}^2-2p_i{p}_j\cos \omega \tilde{\tau }},$$

$$\Im \{S_{ij}\}=\frac{p_i({\mu }_i-{\mu }_j{p}_i{p}_j)\sin \omega {\tau }_i-p_j({\mu }_j-{\mu }_i{p}_i{p}_j)\sin \omega {\tau }_j}{1+(p_i{p}_j{)}^2-2p_i{p}_j\cos \omega \tilde{\tau }}.$$

$$\begin{array}{r}{S}_{ij}(\omega )=\frac{2\lambda p(1+p)\cos \omega \tau }{1+p^4-2p^2\cos 2\omega \tau },\end{array}$$

$$\hat{P}(t,t+\Delta t;t+s,t+s+\Delta t)=\sum _{i=1}^2\sum _{j=1}^2{P}_{ij},$$

$$\hat{C}(s)=\frac{1}{T}{\int }_0^{T}dt\underset{\Delta t\to 0}{\lim }\frac{\hat{P}(t,t+\Delta t;t+s,t+s+\Delta t)}{\Delta t^2}=\sum _{i=1}^2\sum _{j=1}^2{C}_{ij}$$

$$S_X(\omega )=\sum _{i=1}^2{\mu }_i\frac{1-(p_i{p}_j{)}^2+2p_i(\cos \omega {\tau }_i-p_i{p}_j\cos \omega {\tau }_j)}{1+(p_i{p}_j{)}^2-2p_i{p}_j\cos \omega \tilde{\tau }}.$$

$$Q_i(T)=\{\begin{array}{ll}1-e^{-{\mu }_{i}T},& T<\tilde{T},\\ 1-(1-\tilde{P})e^{-{\mu }_i\tilde{T}-{\tilde{\mu }}_i(T-\tilde{T})},& T\ge \tilde{T}.\end{array}$$

$${\tilde{\mu }}_i=$$

$$={\lambda }_i+\sum _{l=1}^{n-1}{\lambda }_{i-l}\prod _{j=1}^l{p}_{i-j}.$$

$${\mu }_i=\underset{m\to \infty }{\lim }{\tilde{\mu }}_i(1+\tilde{P}+{\tilde{P}}^2+\cdots +{\tilde{P}}^m)=\frac{{\tilde{\mu }}_i}{1-\tilde{P}}.$$

Full text

Excitable systems are basically in a resting state, but can generate a strong output under a small but finite forcing. A prominent and very important example in neuroscience is a neuron, which produces a spike when the input forcing is strong enough. Under noisy action, an excitable neuron produces a sequence of random spikes. In this paper we show, that with an additional time-delayed coupling, a network of noisy excitable neurons can produce rather coherent bursts - sequences of spikes separated by fixed delay times; the number of spikes in a burst is random. We construct a point process model for this stochastic bursting, and derive analytically the properties of the interspike intervals, and of the correlations and the spectra.

I Introduction

Processes in delay-coupled nonlinear elements have been attracting a lot of attention, in oscillators (or neurons) Yeung and Strogatz (1999); Rosenblum and Pikovsky (2004); Perlikowski et al. (2010), laser dynamics Soriano et al. (2013), and complex networks D’Huys et al. (2008); Flunkert et al. (2010). While for deterministic oscillators the major interest is in synchronization phenomena, for noise-induced processes time-delayed coupling is known to effect strongly the coherence and the correlation properties of the processes.

There are two basic models of noise-induced processes: noise-induced switchings between two stable states (resulting in a telegraph-type stochastic process), and noise-induced pulses in an excitable system. For the former situation, the cross-correlations between bistable units were investigated in unidirectional delay-coupled networks Kimizuka and Munakata (2009), by extending the theory of a lump bistable noisy model with a delayed feedback Tsimring and Pikovsky (2001). In another approaches, one studies delay-coupled systems with noise within the mean field approximation framework Huber and Tsimring (2003); Hasegawa (2004), focusing on the evolution of the mean and the variance of the order parameter while igonoring the correlation of different units. A full understanding of situations with more complex coupling topology, e.g. for the all-to-all couplingKimizuka et al. (2010), is still an ongoing subject of research.

Noise-induced pulses in an excitable unit is one of the basic models in neuroscience. In this context, delayed feedbacks and couplings are very natural due to a finite time of pulse propagation in connecting synapses. In our previous paper Zheng and Pikovsky (2018), a novel delay-induced spiking pattern which we called stochastic bursting, was observed in a single noisy excitable system with a time-delayed feedback (in the context of neuroscience this corresponds to an autapse with a finite propagation time). This stochastic bursting can be characterized as a random sequence of quasiregular patches, with a pronounced peak in the spectrum at the frequency corresponding to the delay time.

In this paper we extend the theory Zheng and Pikovsky (2018) to the case of two mutually coupled excitable units, and further to a simple, but widely used, topology of an unidirectional delay-coupled network. After outlining the main features and approximations behind the theory of one unit, we describe stochastic bursting in two coupled units in details. The generalization to a chain of neurons will be then straightforward.

II Basic model and properties of one unit

We consider in this paper a network of unidirectionally coupled units, topology of which is illustrated in Fig. 1a. The units are generally different, and the propagation times for the interactions are also different. Each unit is described with a prototypic model for an excitable system, a noisy theta-neuron Ermentrout and Kopell (1986) (or, in other contexts, called active rotator Shinomoto and Kuramoto (1986)):

[TABLE]

Variable $\theta$ is defined on a circle $0\leq\theta<2\pi$. Here parameters $a_{i}$ define the excitability property of the neurons. For $a_{i}<1$, there is a stable and an unstable steady states for an isolated unforced unit, and these states collide in a SNIPER bifurcation at ${a_{i}=1}$. Thus, close to this threshold, the unit is excitable: a small noise or a small force may induce a spike (nearly $2\pi$-rotation of $\theta$ back to the stable state on the circle). Parameter $D_{i}$ describes intensity of the Gaussian white noises $\xi_{i}(t)$, with $\langle\xi_{i}(t)\rangle=0$, $\langle\xi_{i}(t)\xi_{j}(t^{{}^{\prime}})\rangle=2\delta_{ij}\delta(t-t^{{}^{\prime}})$. Finally, parameters $\epsilon_{i}$ describe the strengths of delayed coupling. The coupling force, amplitude of which is $\epsilon_{i}$, is chosen to vanish in the steady state of a driving unit. The forcing term on the r.h.s. of (1) produced by one spike can be represented as

[TABLE]

with $\Theta_{sp}(t)$ being the deterministic trajectory connecting the unstable point (the threshold) with the stable one:

[TABLE]

III Network dynamics and the point process model

To describe qualitatively the dynamics in the network, we start with a non-coupled unit. For a small noise, it produces independent spikes, constituting a Poisson process with rate $\lambda$. Calculation of this rate is a standard task. One formulates the Fokker-Planck equation for the evolution of the probability density of a noisy unit obeying Eq. (1):

[TABLE]

The stationary solution of (4) is

[TABLE]

Here $C$ is a normalization constant. Then the probability current yields the rate of spontaneous spike excitations:

[TABLE]

With coupling, i.e. with $\epsilon\neq 0$, spikes of neuron $i-1$ produce, with a delay, a kick to its next neighbor $i$. Such a kick facilitates excitation, so that it will cause a pulse in neuron $i$, as a combinational effect of forcing and noise, with probability $p_{i}$. The timing of the induced spike is slightly shifted to the forcing, we will denote this shift $\bar{\tau}$ and define a new effective delay $\tau_{i}=\hat{\tau}_{i}+\bar{\tau}$ (we will mainly consider the shift $\bar{\tau}$ as a fixed one, but will briefly discuss the effect of fluctuations of this quantity in section VI.) Hence, to be more general, units are described by different quantities $\lambda_{i}$ (rates with which they “spontaneously” produce spikes due to noise), $p_{i}$ (the probability with which a spike is induced by the incoming force), and $\tau_{i}$ (time delay in forcing).

This allows us to describe the activity on the network as a point process, in which we neglect the durations of the spikes (approximate them as delta-functions), compared to the delay times $\tau$ and the inverse rate $\lambda^{-1}$. This is well justified for mammal brains, where the characteristic duration of a spike is $\sim 1ms$, while the delay time and the characteristic time interval between noise-induced spikes are of order $\sim 100ms$ Swadlow and Waxman (2012); Longtin (2013). The spike occurred at time $t$ in neuron $i$, will produce a kick on neuron $i+1$ at time $t+\hat{\tau}_{i}$, and will generate a spike in neuron $i+1$ at time $t+\tau_{i}$ with probability $p_{i}$, as depicted schematically in Fig. 1 (a).

Throughout the paper, in numerical illustrations we use parameters $a=0.95$, $D=0.005$. The small additional delay is $\bar{\tau}\approx 7$, it is much smaller than the characteristic delay times we use $\hat{\tau}\gtrsim 100$ and the inverse of the spontaneous rate $\lambda^{-1}\approx 1506$. For these parameters of the neuron, we use the coupling strength $\epsilon=0.14$, for which the probability to induce a spike by the forcing is $p=0.53$ (for details of the calculation of this probability see Ref. Zheng and Pikovsky, 2018). Empirically, this probability can be determined in simulations of one unit with a delayed self-feedback. One calculates the numbers of spikes, during a large time interval, in dependence on the coupling strength $N(\epsilon)$. Then,

[TABLE]

IV Two coupled units

IV.1 Statistics of interspike intervals

It is instructive to start with the case where there are only two neurons in the ring, i.e. $n=2$, and then to extend the theory to a more general case with arbitrary $n>2$. When $n=2$, the two delay-coupled neurons are denoted as $i$ and $j$ ($i,j$ could be $1$ or $2$). We simulate Eq. (1) and obtain spike trains with bursts of each neuron as shown in Fig. 1(b). As outlined above, an idealized point process can be constructed to describe the bursting phenomenon, as illustrated in Fig. 2. Spontaneously generated spikes we denote as leaders. Each leader induces a finite set of followers (induced spikes), and together with them constitutes a coherent burst. In a burst, spikes in the two neurons appear alternatively with time intervals $\tau_{1}$ and $\tau_{2}$. The number of spikes in a burst is random. Noteworthy, similar to the case of one neuron with delayed feedback Zheng and Pikovsky (2018), the bursts can overlap; hence the analysis of the ISI (inter-spike intervals) distribution of the spike train in each unit is nontrivial. Compared to the single neuron case Zheng and Pikovsky (2018), here the leaders in each neuron will have random followers in both neurons, as explained schematically in Fig. 2.

First, we determine the overall rate of the spikes in each unit. The probability for a spike in unit $i$ to induce a follower in the same unit is $p_{i}p_{j}$. Thus, the probability for a leader to have exactly $m$ followers in the same unit is $(p_{i}p_{j})^{m}(1-p_{i}p_{j})$. The average number of followers in the same unit is $\sum_{m}m(p_{i}p_{j})^{m}(1-p_{i}p_{j})=(p_{i}p_{j})(1-p_{i}p_{j})^{-1}$. The total average number of spikes in a burst is $1+(p_{i}p_{j})(1-p_{i}p_{j})^{-1}=(1-p_{i}p_{j})^{-1}$ Therefore, the total rate of spikes initiated in unit $i$ is $\lambda_{i}(1-p_{i}p_{j})^{-1}$. For the spikes in unit $i$, initiated by a leader in unit $j$, we have first to find the rate of the first followers in unit $i$, which is $\lambda_{j}p_{j}$; the total rate of these spikes is thus $\lambda_{j}p_{j}(1-p_{i}p_{j})^{-1}$. Summing, we obtain the total rate of spikes $\mu_{i}$ as

[TABLE]

To derive the statistics of the ISI, we assume that in one of the units there is a spike at time $t$, and the next spike at time $t^{\prime}>t$, so that the inter-spike interval is $T=t^{\prime}-t$. Three different cases should be distinguished, namely, $T>\tau_{i}+\tau_{j},T=\tau_{i}+\tau_{j}$, and $T<\tau_{i}+\tau_{j}$. If $T<\tau_{i}+\tau_{j}$, the spikes at $t$ and $t^{\prime}$ can be either spontaneous (leader) or delay-induced ones, but in the latter case they belong to different bursts, so they are independent. Therefore, the survival function, i.e., the probability that there is no spike in $(t,t^{\prime})$, is determined by the full rate $\mu_{i}$ from (8): $S_{i}(T)=\exp(-\mu_{i}T)$.

In contradistinction, for the case $T>\tau_{i}+\tau_{j}$, the spike at time $t^{\prime}$ in neuron $i$ can only be from a spontaneous one (leader) in neuron $i$ itself, or the first induced spike (with probability $p_{j}$ from neuron $j$). These events are independent on the occurrence of a spike at time $t$, and have the total rate $\lambda_{i}+\lambda_{j}p_{j}$. The probability that there is no any spike in $(t,t^{\prime})$ in neuron $i$ is the product of three terms: the probability to have no spikes in the interval $(t,t+\tau_{i}+\tau_{j}]$ with the survival function $S_{\tau b}=\exp(-\mu_{i}(\tau_{i}+\tau_{j}))$, the probability $(1-p_{i}p_{j})$ not to have a follower for the spike at $t$, and the probability to have no spike in the interval $(t+\tau_{i}+\tau_{j},t^{\prime})$, where only the spontaneous total rate $\lambda_{i}+\lambda_{j}p_{j}$ applies with the survival function $S_{\tau a}=\exp(-(\lambda_{i}+\lambda_{j}p_{j})(T-\tau_{i}-\tau_{j}))$. Thus, the survival function for the case $T>\tau_{i}+\tau_{j}$ is

[TABLE]

Based on the above description, and on the relationship between the cumulative ISI distribution $Q(T)$ and the survival function $S(T)$, which reads $Q(T)=1-S(T)$, the cumulative ISI distribution of neuron $i$ can be obtained as follows:

[TABLE]

where $\tilde{\tau}=\tau_{i}+\tau_{j}$. We compare this expression with results of numerical simulations in Fig. 3. The point process described by Eq. (10) agrees well with direct simulations of Eq. (1), where we assumed the parameters of the two neurons to be the same, except for the time delays which are different.

IV.2 Correlations and spectra

In the following we derive the autocorrelation and the cross-correlation functions of the spike trains, and the corresponding power spectrum and the cross-spectrum. The autocorrelation function is defined via a joint probability to have a spike in unit $i$ within a small time interval $(t,t+\Delta t)$, and a spike in unit $j$ within the time interval $(t+s,t+s+\Delta t)$, no matter whether or not there are any spikes between $t$ and $t+s$. The joint probability of these events is defined as

[TABLE]

Here $W_{i}(t,t+\Delta t)=\mu_{i}\Delta t$ is the probability to observe a spike in neuron $i$ within the time interval $[t,t+\Delta t]$. The quantity $P_{ij}(t+s|t,\Delta t)$ is the probability to induce a spike in neuron $j$ at time $t+s$, given a spike in neuron $i$ at time $t$.

IV.2.1 Correlations and spectra within one unit

We first calculate the conditional probability (11) for the same unit. The conditional probability to have one induced spike is $p_{i}p_{j}$, this event happens with time shift $\tilde{\tau}$; the conditional probability for the $k$-th induced spike is $(p_{i}p_{j})^{k}$, it happens with delay $k\tilde{\tau}$. Therefore,

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function. Since the correlation function can be seen as the mean rate of the joint event, after substituting Eq. (11) and (12), the auto-correlation function is

[TABLE]

where we use the fact that $P_{ii}$ is $t$-independent, and have taken into account that $\langle x_{i}\rangle=\mu_{i}$.

Taking the Fourier transform of the correlation function, we obtain the power spectral density

[TABLE]

The derivations above are based on the time series $x_{i}(t)$ represented as a sum of delta-peaks, i.e., $x_{i}(t)=\sum\limits_{i=1}^{N}\delta(t-t_{i})$. For a train of realistic spikes, the shape function can be straightforwardly taken into account as done in Ref. Zheng and Pikovsky, 2018, namely the spectrum (14) should be just multiplied by the squared amplitude of the Fourier transform of the pulse shape. For example, if observable (2) is used, the spike train $x_{i}$ will be convoluted with the shape function $H(t)$. Hence, the power spectral density and the cross-spectral density in the following illustrations will be multiplied by the spectral density of $H(t)$, which we denote as $S_{H}(\omega)$. For simplicity, in the formulas below we still use the delta-peak representation of $x_{t}$, while we multiply by $S_{H}(\omega)$ to compare with numerical correlations and spectra obtained by simulations of Eq. (1). This comparison is shown in Fig. 4(a), the theoretical predictions based on the point process analysis agree well with the results of direct simulation of Eq. (1).

IV.2.2 Cross-correlations and cross-spectra for two units

The conditional probability of the joint event between the two neurons can be expressed similarly to formula (12) above:

[TABLE]

This allows us to obtain the cross-correlation function of the two neurons, by substituting Eq. (15) into Eq. (11), leading to the Eq. (16):

[TABLE]

The cross-spectral density is the Fourier transform of the cross-correlation function:

[TABLE]

It is instructive to present explicitly the real part

[TABLE]

and the imaginary part

[TABLE]

of the cross-spectrum.

Unlike the power spectral density described by a real-valued function (14), the cross-spectral density is generally a complex-valued function. It is real-valued only when the two neurons are totally identical, i.e., $\lambda_{i}=\lambda_{j}=\lambda$, $p_{i}=p_{j}=p$, and $\tau_{1}=\tau_{2}=\tau$, resulting in a simple expression

[TABLE]

which is very similar to the power spectral density of a single unit (14). We compare the theoretical cross-spectra with simulations in Figs. 4b,c.

IV.2.3 Correlation and spectra of the total output from the network

If we consider correlations and spectra from the viewpoint of the total network output, the cross-correlations between all the pulses should be calculated. A joint probability could be defined as having a spike in any unit within a small time interval $(t,t+\Delta t)$, and a spike in any unit within the time interval $(t+s,t+s+\Delta t)$, no matter whether or not there are any spikes between $t$ and $t+s$. The joint probability $\hat{P}$ of these events is the sum of all contributions:

[TABLE]

where $P_{ij}$ is described by Eq. (11). Thus the correlation function is

[TABLE]

where $C_{ij}$ is described by Eq. (13) when $i=j$ and by Eq. (16) when $i\neq j$. The spectral density of the total output , i.e. of the observable $X(t)=x_{1}(t)+x_{2}(t)$, is obtained as a Fourier transform of $\hat{C}(s)$, leading to

[TABLE]

As is shown in Fig. 5 (a), the theoretical spectra of the total output agrees well with simulation results.

V General network

The case of many neurons with $n>2$ in the ring topology is a direct extension of the $n=2$ case as described above; thus the analysis follows the same steps, only the expressions are more involved. First, we extend the cumulative ISI distribution for neuron $i$ in the ring as follows,

[TABLE]

Here $\tilde{T}=\sum\limits_{i=1}^{n}\tau_{i}$ is the total round-trip delay time acround the ring, $\tilde{P}=\prod\limits_{j=1}^{n}p_{j}$ is the probability to have a completed round trip around the ring, and $\tilde{\mu}_{i}$ is the spiking rate of all first spikes in bursts that include neuron $i$:

[TABLE]

Here $\lambda_{i}$ is the rate of spontaneous spikes in neuron $i$ itself, $\lambda_{i-1}p_{i-1}$ is the rate of spontaneous spikes in neuron $i-1$ that induce also a spike in neuron $i$, and so on. Noteworthy, due to the ring structure of the coupling, $\lambda_{i}$ is a periodic function, i.e., $\lambda_{i}=\lambda_{i+n}=\lambda_{i-n}$. The total activity in the network is characterized by the rate $\mu_{i}$, to which the rates from all spikes (both leaders and followers) contribute. Thus, similarly to Eq. (8), we get

[TABLE]

Differently formulated, the expression above follows from the fact that a spike can have $m$ followers (in the same unit) with probability $\tilde{P}^{m}(1-\tilde{P})$.

Using the same method as in the $n=2$ case described above, we obtain the power spectral density of neuron $i$ in the ring by Fourier transform of the correlation function (not presented):

[TABLE]

The cross-spectral density of spike trains in neuron $i$ and neuron $j$ is:

[TABLE]

The real part of which is:

[TABLE]

and the imaginary part of which is

[TABLE]

Here $T_{ij}=\tau_{i}+\cdots+\tau_{j-1}$ is the delay time from neuron $i$ to neuron $j$ along the direction of the ring , i.e. clockwise as depicted in Fig. 1, with probability $\bar{P}_{ij}=\prod\limits_{l=i}^{j-1}p_{l}$. Correspondingly, $T_{ji}=\tilde{T}-T_{ij}$ is the delay time from neuron $j$ to come to neuron $i$ with probability $\bar{P}_{ji}$ and $\bar{P}_{ij}\bar{P}_{ji}=\tilde{P}$. The spectral density of total output, i.e. $X=\sum\limits_{i=1}^{n}x_{i}(t)$, from the network is

[TABLE]

In the case that all the units are totally identical, i.e. $\lambda_{i}=\lambda$, $p_{i}=p$ and $\tau_{i}=\tau(i=1,\cdots,n)$, $S_{X}(\omega)$ reduces to $S_{X}(\omega)=\frac{n\lambda(1+p)}{1+p^{2}-2\cos\omega\tau}$.

Generally, the model works for any network size $n$, but for simplicity we choose $n=3$ for comparison with numerics. The cumulative ISI described by Eq. (24), spectra described by Eqs. (27), (29), (30) and (31) agree well with direct simulation of Eq. (1), as shown in Fig. 3(b), Fig. 4(d)-(f) and Fig. 5(b), respectively. Noteworthy, similar to the case $n=2$, the cross-spectrum $S_{ij}$ is generally a real-valued function only if $n$ is an even number and neurons $i$ and $j$ are symmetric, i.e., $|i-j|=n/2$.

To further demestrate that our theory works for a larger network, we choose $n=10$ and calculate the cross-spectra between neurons at different distances, e.g. between neurons 1 and 3, and between neurons 1 and 4. As shown in Fig. 6, the analytical results agree well with the simulations. Noteworthy, as $\epsilon$ goes larger, the duration of the delay-induced pulse becomes shorter, leading to a smaller empirical time shift $\bar{\tau}$. In the case depicted in Fig. 6, $\bar{\tau}\approx 5$ for $\epsilon=0.2$.

VI Conclusion

In conclusion, we investigated the stochastic bursting phenomenon in $n$ unidirectional delay-coupled noisy excitable systems. Under the condition of time-scale separation, an idealized version of coupled point processes with leader-follower relationship was formulated. Roughly speaking, occurrence of stochastic bursting is based on three ingredients: excitability of the system, excitatory coupling with a fixed time delay, and noise. Excitability combined with noise results in the spontaneous spikes with a constant spiking rate, which are leaders of the bursts. A relatively weak coupling is not strong enough to induce a follower deterministically, but it leads to an increased probability to have a follower, characterized by the crucial parameter $p$. The leader with the followers form a burst, which is rather coherent (because of the fixed time interval between the followers, nearly equal to the delay time), but has a random number of spikes in it.

To characterize the stochastic bursting, the cumulative ISI distribution was derived; simulations demonstrated a good agreement with the theoretical prediction. Furthermore, via the calculation the joint probability of the spikes, both the auto-correlation function of a single neuron spike train, the cross-correlation function of any pair of neurons in the unidirectional ring, and the auto-correlation function of the total output from the network are derived analytically. Calculation of the spectra and of the cross-spectra is then straightforward. Noteworthy, the model in the present paper not only shows an interesting coherent spiking pattern, but also provides an alternative way to investigate the cross-spectrum of different neurons beyond the linear response theory (see, e.g., Refs. Lindner et al., 2005; Ostojic et al., 2009; Trousdale et al., 2012; Vüllings et al., 2014, to name a few), which is widely used in the analysis of correlated neuronal networks.

Above we assumed, based on the time scale separation, that the delay times are constants. A generalisation to the case of random delay times is also possible and will be presented in details elsewhere; here we discuss a simple version of this analysis. The essential point where the fixed delays appear, is the representation of the correlation function (13) as a sum of delta-peaks at times $n\tilde{\tau}$. If we assume the delay times to be independent Gaussian variables with mean value $\tilde{\tau}$ and standard deviation $\kappa$, then one has to replace in (13) delta-functions by Gaussian peaks $\delta(t-n\tilde{\tau})\to(2n\pi\kappa^{2})^{-1/2}\exp[-(t-n\tilde{\tau})^{2}/(2n\kappa^{2})]$. In the spectrum  (14), this correction corresponds to the replacement $p_{i}p_{j}\to p_{i}p_{j}\exp[-n\omega^{2}\kappa^{2}/2]$. Around the main frequency peaks (i.e. with small values of $n$), the effect of this correction is, as expected, small, due to the time scale separation $\kappa\ll\tilde{\tau}$.

In this paper we restricted our attention to a unidirectional coupling in the ring geometry, because here overlapping of incoming spikes is not possible (or, better to say, is very unprobable under the condition of the time-scale separation). Such an overlap happens, e.g., in a network of delay-coupling neurons demonstrating polychronization Izhikevich (2006); study of stochastic bursting in such a setup is a subject of ongoing research.

Acknowledgements.

C.Z. acknowledges the financial support from China Scholarship Council (CSC). A.P. was supported by the Russian Science Foundation (Grant No. 17-12-01534).