Stochastic bursting in networks of excitable units with delayed coupling

We investigate the phenomenon of stochastic bursting in a noisy excitable unit with multiple weak delay feedbacks, by virtue of a directed tree lattice model. We find statistical properties of the appearing sequence of spikes and expressions for the power spectral density. This simple model is extended to a network of three units with delayed coupling of a star type. We find the power spectral density of each unit and the cross-spectral density between any two units. The basic assumptions behind the analytical approach are the separation of timescales, allowing for a description of the spike train as a point process, and weakness of coupling, allowing for a representation of the action of overlapped spikes via the sum of the one-spike excitation probabilities.


Introduction
Bursting, which plays an important role in neuronal communication and synchronization, refers to a dynamical state where a neuron repeatedly fires a relatively regular finite sequence spikes; bursts are separated by epochs where the neuron is in a resting state (Izhikevich 2000;Coombes and Bressloff 2005). Bursting is observed in neurons of different types, such as in neocortex (Connors and Gutnick 1990), hippocampus (Dzhala and Staley 2004;Su et al. 2001) and cerebellum (Womack and Khodakhah 2002); see the classification in (Izhikevich 2006b). In many situations, bursting is an intrinsic property of a neuron, following from the particular properties of its native dynamics. Correspondingly, the theory of such bursting is based on the bifurcation theory of Communicated by Peter Thomas.
In our previous works Pikovsky 2018, 2019), we have demonstrated that a coherent spike pattern, which we call stochastic bursting, can appear in simple excitable units due to the combined effect of time-delayed feedback and noise. Noise plays a twofold role in the stochastic bursting. On the one hand, it leads to spontaneous appearance of quite rare statistically independent spikes. On the other hand, when a delayed feedback pulse enters like an excitatory underthreshold kick (i.e., a weak kick that itself without noise does not produce a new spike), noise results in an increased probability to create a new induced spike. Thus, a spontaneous spike (the leader) can be followed by a sequence of induced spikes (the followers) separated approximately by the delay time interval. Because the creation of a follower is a random event due to noise, the overall burst has a random number of spikes. We described stochastic bursting statistically in the case of a single excitable unit (Zheng and Pikovsky 2018) and for networks of unidirectionally delay-coupled units (Zheng and Pikovsky 2019). What these two cases have in common, is that any two delay-induced kicks do not overlap. This allowed for a full statistical description of the bursting as a point process, where the only parameters are the spontaneous rate of excitation and the probability to excite a follower. The point process model is an idealization based on the timescale separation: It is assumed that the characteristic duration of Fig. 1 Schematic description of a noisy excitable unit with two feedback loops a spike is much less than other characteristic times in the system, the delay time and the characteristic time interval between the pulses (which depends on the spontaneous rate and the probability of induced spikes).
However, it becomes more challenging when neurons have multiple feedback or more complex coupling topology, where two or more delay-induced spikes could overlap. Such an overlap of incoming spikes in a purely deterministic setup let E. Izhikevich to introduce the concept of polychronization (Izhikevich 2006a). Specific polychronous patterns are determined by combinations of delay time. They appear if, for a reliable excitation, superposition of delay-induced kicks is favorable. Here, below we treat the effect of overlap of incoming kicks in a stochastic setting, assuming that the probability to induce a spike increases if two kicks overlap.
Similar to the previously studied cases of stochastic bursting (and generally for systems with noise and delay), the process under study is non-Markovian, which prevents using the powerful Fokker-Planck formalism to obtain analytic results. Thus, we adopt the point-process representation as a good option to describe the statistics of the spike train in each neuron and between different neurons, if timescale separation holds.
In the present paper, we extend the theory of stochastic bursting to the case where incoming delayed pulses can overlap. Here, generally, one needs to define the probability to induce a spike by two incoming pulses. We restrict our attention below to the case of weak coupling, where this probability can be represented through one-pulse probabilities. This allows for an analytical treatment, which results in explicit expressions for the power spectrum of the spike trains. First, in Sect. 2, we investigate the stochastic bursting phenomenon in a noisy excitable unit with multiple delay feedbacks. Then, in Sect. 3, we extend the theory to a network of three delay-coupled units, with a star-type coupling.

One excitable unit with multiple delayed feedbacks
The main goal of this paper is to extend the theory of stochastic bursting to the situations where an overlap of delayed feedback pulses is possible. The simplest case is one unit with two delayed feedbacks, as shown in Fig. 1. (In one unit with one delay, the incoming pulses cannot overlap.) As we see in Sect. 3, the hub unit (unit 2 in Fig. 6) in a star-type network possesses two or several feedback loops with different delays and is thus similar to the simplest setup of this section. As a model, we consider a scalar equation on a circle 0 ≤ θ < 2π , which is a prototypical example of excitability: (1) Here, parameter a 1 defines the excitability level and parameters 1 and 2 are coupling strengths for the two delayed feedbacks with delay timesτ 1 andτ 2 , respectively. The system is driven by a Gaussian white noise ξ(t) with intensity D, i.e., ξ(t) = 0, ξ(t)ξ(t ) = 2Dδ(t − t ). By direct simulation of the Langevin equation (1) using Euler-Maruyama scheme with time step Δt = 0.01, we observe rather coherent spiking which we call stochastic bursting in the spike train, as shown in Fig. 2a.
Without feedback (i.e., for 1 = 2 = 0), the model is equivalent to the well-known models of active rotators (Shinomoto and Kuramoto 1986) and to the popular theta model (Ermentrout and Kopell 1986). For a < 1, there are two steady states, one stable and one unstable. Noise drives the state of the system beyond the unstable steady state, and a rotation around the circle back to the stable state occurs; this rotation is indicated as a spike. To find statistics of the spikes, one formulates the Fokker-Planck equation for the evolution of the probability density of a noisy unit obeying Eq. (1) as follows: The stationary solution of (2) iŝ Here, C is a normalization constant. The probability current across threshold yields the rate of spontaneous spike excitations: Below we will assume this rate to be small, i.e., the characteristic time interval between the noise-excited spikes is much larger than the duration of the spikes. Below we will also assume that this separation of times is valid for characteristic times of the delayed feedback: The delay timesτ 1 ,τ 2 and their difference |τ 1 −τ 2 | are much larger than the duration of a pulse.  (1) with two feedback signals, with values of parameters chosen as a = 0.95, D = 0.005, 1 = 0.12,τ 1 = 500, 2 = 0.1 andτ 2 = 600. Zooming in a burst shows many spikes with interspike interval of around 500 and 600. (b) Representation in directed tree lattice, where the cross-position of the lattice, i.e., kτ 1 +lτ 2 with k, l being non-negative integers, represents where there is potentially a delay-induced spike with some probability P(kτ 1 + lτ 2 ). The blue numbers near the intersection kτ 1 + lτ 2 are the number of paths to go to that point, i.e., k+l k in Eq. (7) We now qualitatively describe the phenomenon of stochastic bursting illustrated in Fig. 2. If just one delay feedback term is present, the bursting appears as described in Zheng and Pikovsky (2018). Because, after a pulse, a feedback force acts as a kick on the unit after delay timeτ , there is an increased probability p for a next pulse to be induced. Thus, each spontaneously excited pulse can have a random number of "followers" separated by a time interval τ =τ + τ r ; together, they constitute a regular burst. Here, we take into account that the delay-induced effect is not instant, but rather it takes some relatively small response time τ r for the unit to generate a spike after receiving a delayed kick, as described in Ref. Zheng and Pikovsky (2019). Generally speaking, the delay-induced spike occurs with some deviation around the effective time delay, i.e., not exactly at t +τ + τ r , but for simplicity, we assume the deviation can be ignored. The essential parameter of the stochastic bursting is the probability p to induce a follower, and we discussed in Ref. Zheng and Pikovsky (2018) the ways to calculate it.
In the case of two delayed feedback signals, after a spontaneous pulse at time t, there are two times t + τ 1 and t + τ 2 at which the followers can appear independently, with probabilities p 1 and p 2 , and these probabilities are the same as for single time-delay feedbacks, as shown schematically in Fig. 2b. However, at the next level, there is an interaction of delayed kicks, which makes the problem essentially more difficult than that of one delayed feedback. Below in this paper we will assume that times τ 1 and τ 2 are incommensurate. Thus, we exclude resonances like τ 1 = 2τ 2 , which will be treated in a separate work. We stress here that practically, because we consider relatively weak feedback, the number of spikes in a burst is not large. This means that "incommensurability" should be understood in a weak sense-as absence of resonances mτ 1 = nτ 2 with small integers m, n.
If both followers at times t + τ 1 and t + τ 2 are present, at time t + τ 1 + τ 2 there will be an overlap of two incoming feedback kicks. Generally, there are many such overlaps possible, at times t + kτ 1 + lτ 2 with k and l being both positive integers. (Here, due to incommensurability mentioned, for each overlap there is only one pair of the integers k, l which yields it.) Thus, another probability p (2) to induce a spike mediated by the overlap of two incoming kicks is needed. Generally, one must calculate this probability de novo. One may employ methods similar to those used for determining p 1 and p 2 . However, when the two feedbacks are both weak and independent, we can assume a linear dependence of the probability to induce a follower on the incoming pulse amplitude, which makes an approximation p (2) ≈ p 1 + p 2 , to be adopted below, reasonable.
It is easy to see that potential followers of a spontaneous spike form a tree, as illustrated in Fig. 2b. On this tree, there can be non-branching paths, which correspond just to sequences of followers separated by time intervals τ 1,2 , which appear with probabilities p 1,2 . However, any branching leads to an overlap, so we use p (2) to calculate the probability to observe the corresponding induced spike. We note here that the considered stochastic model is not a standard branching process, because of the overlaps.
It is instructive to start with the simplest overlap at time τ 1 +τ 2 . The probability to induce a spike at t = τ 1 +τ 2 given a spike at t = 0, i.e.,P(τ 1 + τ 2 ), can be calculated as

Fig. 3
Configurations and the corresponding probability to induce a spike at the cross section τ 1 + τ 2 (a)-(c) and 2τ 1 + τ 2 (d)-(j) of the lattice shown in Fig. 2b Here, the first line in Eq. (5) describes contributions of different configurations with corresponding probabilities, as shown in Fig. 3a-c. Similarly, we can calculate the probability to induce a spike at the time-delayed instant 2τ 1 + τ 2 as The corresponding configurations and their probabilities are shown in Fig. 3d-j. Now, by induction, it is easy to extend to the general case. One can easily check that the general expressioñ is consistent with calculation of the induced probabilitỹ P(kτ 1 + lτ 2 ) on the base of the "parent" probabilities P((k − 1)τ 1 + lτ 2 ) andP(kτ 1 + (l − 1)τ 2 ). By iteration of the relationship on the directed tree lattice, it is easy to obtain Eq. (7). We stress that expression (7) Fig. 4 Spike rate of the single excitable unit with two delayed feedbacks, where 1 is set to be 0.12 and 2 is varied. The blue circles represent the result by direct simulation of Eq. (1), and the red dashed line represents the analytical result described by Eq. (9). Here, τ 1 and τ 2 are chosen as 500 and 600, respectively p (2) ≈ p 1 + p 2 ; in a general case p (2) = p 1 + p 2 , we could not derive a simple expression for these probabilities.
Having determined the probabilities of the followers, we now derive statistical properties of the bursts. To calculate the total firing rate μ, we have to determine the average number of spikes in a burst. The probability to have k spikes with τ 1 separation and l spikes with τ 2 separation in the burst is the product ofP(kτ 1 +lτ 2 ) and probability 1− p 1 − p 2 , which is the probability to generate no spikes further. Thus, the total firing rate is We check this relation numerically in Fig. 4. Throughout the paper, the values of parameters for a and D are fixed as a = 0.95 and D = 0.005. This yields the spontaneous spike rate λ = 6.64 × 10 −4 , as calculated from Eq. (4). The probabilities to induce a spike can be calculated by virtue of solving a forced Fokker-Planck equation numerically as described in Ref. Zheng and Pikovsky (2018). For the delay coupling amplitudes = 0.1 and = 0.12, this gives values p = 0.25 and p = 0.39, respectively. Furthermore, in this case, the empirical value of the response time τ r is approximately τ r ≈ 7. In Fig. 4, we set the value of 1 fixed, i.e., Our next goal is to calculate correlations and spectra of stochastic bursting in this system. The autocorrelation function C(s) can be represented in terms of the joint probability to have a spike in the time interval (t, t + Δt) and a spike in the time interval (t + s, t + s + Δt): Since the joint probability P(t, t +Δt; t +s, t +s +Δt) is the product of the probability to fire a spike in the time interval [t, t + Δt], i.e., μΔt , and the conditional probability to fire a spike within the time interval [t + s, t + s + Δt] given a spike at time t (see also Ref. Zheng and Pikovsky 2019), i.e., k,lP (kτ 1 + lτ 2 )δ(s ± (kτ 1 + lτ 2 ))Δt, we obtain by using expression (7) The Fourier transform of the correlation function gives the power spectrum: We compare this theoretical prediction with the results of numerical simulation in Fig. 5a. Using the same method, a noisy excitable unit with m delays described by the following Langevin equatioṅ θ = a + cos θ + 1 (a + cos θ(t −τ 1 )) + · · · + m (a + cos θ(t −τ m )) + ξ(t) can be studied. Here, the total spike rate can be expressed as (1) (blue lines), the parameters are chosen as 1 = 0.12,τ 1 = 500, 2 = 0.1,τ 2 = 600 in two-delay case and 1 = 2 = 3 = 0.1,τ 1 = 300,τ 2 = 430,τ 3 = 500 for Eq. (13) in three-delay case. The red lines are the analytic results described by Eq. (12) for two-delay case and Eq. (13) for three-delay case. The power spectral density of the spike train is multiplied by the power spectral density of the shape function (color figure online) Fig. 6 Schematic description of three delay-coupled noisy excitable units in a chain, where p i j is the probability to induce a spike due to delay feedback with strength i j and time delay τ i j formula: where τ l =τ l + τ r . Since the condition of timescale separation becomes harder to fulfill for large m (one needs the time intervals between pulses to be large), practically only the case m = 3 is tested in Fig. 5b. As shown in Fig. 5a, b, the analytic results described by Eq. (12) for two delayed feedbacks and Eq. (14) for three delayed feedbacks agree well with the direct simulation of the Langevin Eqs. (1) and (13), respectively. central unit 2 (a hub) coupled to peripheral units 1 and 3. As is clear from Fig. 6, there are four delay times and four probabilities to induce a follower.
The Langevin equations, describing this network, reaḋ θ 1 = a + cos θ 1 + 21 (a + cos θ 2 (t −τ 21 )) + ξ 1 (t), θ 2 = a + cos θ 2 + 12 (a + cos θ 1 (t −τ 12 )) Here, i j (i, j = 1, 2, 3) is the delay feedback strength from unit i to unit j with delay timeτ i j , and In the absence of delay feedback, i.e., i j = 0, the three units fire spontaneously with constant rates λ i (i = 1, 2, 3), described by Eq. (4). For simplicity, the noise intensities here in the three units are chosen to be the same. When the delayed feedback is included, i.e., i j = 0, each spontaneous spike in unit 1 will induce another spike in unit 2 with probability p 12 after time delay τ 12 , and the induced spike in unit 2 will generate spikes in unit 1 with probability p 21 after time delay τ 21 and in unit 3 with probability p 23 after time delay τ 23 . Here, τ i j includes the response time, i.e., τ i j =τ i j + τ r .
Thus, each spontaneous spike, which plays a role of leader, is followed by random number of induced spikes (followers) across the network. Noteworthy, the above description is based on timescale separation as described in Sect. 2. The relation to the model of one excitable unit with two feedbacks described in Section 2 is evident when one considers effective feedbacks from unit 2 to itself: There are two effective delay feedback channels, one with the probability p 21 p 12 and time delay τ 21 + τ 12 to induce a spike back into unit 2 itself through unit 1 and the other one through unit 3 with probability p 23 p 32 and time delay τ 23 + τ 32 . Therefore, the total spike rate of unit 2 can be represented as Here, the numerator λ 2 + p 12 λ 1 + p 32 λ 3 represents the total rate of first spikes in a burst in unit 2. These are spikes initiated spontaneously in unit 2 (rate λ 2 ), and the first followers in unit 2 of spikes spontaneously created in units 1 and 3 (rates p 12 λ 1 and p 32 λ 3 , respectively). The denominator in (16) comes from the same summation as in Eq. (9). Similarly, we can express the total spike rate μ 1 of neuron 1 as μ 1 = λ 1 + μ 2 p 21 = λ 1 + λ 2 + p 12 λ 1 + p 32 λ 3 1 − ( p 21 p 12 + p 23 p 32 ) since the leading spikes in unit 1 are either the spontaneous ones created in unit 1, or induced from total spikes of unit 2. For unit 3, the total rate is Using the same method as described in Sect. 2, we can write the power spectral density of a spike train in unit 2 as wherep 1 = p 21 p 12 ,τ 1 = τ 21 + τ 12 ,p 2 = p 23 p 32 and τ 2 = τ 23 + τ 32 . Similarly, taking into account all the delayed feedback loops connecting neuron 1, we obtain the power spectral density of neuron 1 as follows: Here, we denote T =τ 1 +τ 2 in the first line, and the term p 1p2 e iωT n>0 (p 2 e iωτ 2 ) n is due to the summation of all the feedback loops starting from neuron 1, connecting neuron 2 and 3 and then going back to neuron 1. By the same method, the power spectral density of neuron 3 is We compare these expressions with the results of direct numerical simulations of model (15) in Fig. 7a-c. For networks of coupled units, it is instructive to calculate cross-correlations and cross-spectra. The cross-correlation between neuron 1 and 2 can be described in terms of the joint probability P 12 (t, t + Δt; t + s, t + s + Δt) that there are a spike in the time interval (t + Δt) in neuron 1 and a spike in the time interval (t + s, t + s + Δt) in neuron 2 (see also Ref. Zheng and Pikovsky 2019), i.e., Here,P(kτ 1 +lτ 2 ) is the probability to induce a spike at time kτ 1 + lτ 2 after a spike, and according to Eq. (7), it is Substituting Eq. (23) into Eq. (22), we obtain the crosscorrelation function The cross-spectral density between neurons 1 and 2 is just Fourier transform of the cross-correlation function, i.e., As shown in Fig. 8 (see panels (a) and (b) for the real and the imaginary parts of the cross-spectral density), the analytic results agree well with direct simulation of Eq. (14).  Fig. 8 Real part (a) and imaginary part (b) of the cross-spectral density between neurons 1 and 2 in a chain of three delay-coupled noisy excitable neurons. Values of parameters are chosen the same as in Fig. 7. Noteworthy, the cross-spectral densities are also multiplied by the power spectral density of the shape function S H as in Fig. 7. Blue lines: theory; red lines: analytic results (color figure online)

Conclusion
In summary, we investigated stochastic bursting in a single noisy excitable unit with multiple feedbacks and in a startype network of three delay-coupled units. Both systems are the simplest examples of a neural network with an overlap of incoming spikes. In the deterministic case, this leads to polychronization (Izhikevich 2006a); in this sense, stochastic bursting in such a network may be referred to as noisy polychronization.
Our analysis is based on two approximations. One is that of timescale separation, valid for weak noise and large delay times. It allows us for modeling the process as a point one, so that only time instants of spikes are relevant for correlations and spectra. Another approximation is based on the assumption that the induced probabilities are small, and the probability for two overlapping inputs to induce a spike can be represented as a sum of the corresponding one-input probabilities. This latter assumption appeared to be extremely helpful, as it allowed us to express the probabilities of induced spikes in a simple closed form. We would like to point out that our analysis does not rely specifically on the theta model but is generally applicable to any excitable unit as long as the assumptions above (most important is that of separation of timescales) are met. Once the spontaneous firing rate and probabilities of generating of followers are known, the calculations of the total firing rate of each unit and of the pairwise correlations are straightforward.
As a result of our analysis, the spectra of the point stochastic bursting process have been analytically represented in a closed form. These power spectral and cross-spectral densities (12), (14) and (20-21) agree well with direct simulation of the original Langevin equations (1), (13) and (15). A general structure of the analytic expressions is the same: It yields a power spectrum with peaks at the delay times; the strengths of these peaks are roughly proportional to the corresponding probabilities to create a noise-induced follower. Remarkably, the cross-spectrum (25) has not a peak form, but rather looks like modulated oscillations, where characteristic modulation periods are the partial delays between the units.
Our approach works well for small numbers of interacting units and small numbers of feedback loops, because here the timescale separation is well justified. In a network with a large number of units and many connections, the spikes become denser, and the point process approximation might be violated. Similarly, for strong feedbacks the probability to induce a spike by overlapping inputs will be a non-trivial quantity (not a sum of one-input probabilities) as well.
At this point, we would like to mention several relevant studies of correlations in noisy neural networks. Pernice et al. (2011) investigated the pairwise correlation in networks of neurons with both excitatory and inhibitory interaction, based on the linear Hawkes processes. Trousdale et al. (2012) analyzed how the network structure influences the correlations based on the linear response theory (Lindner et al. 2005). There are also more sophisticated methods to investigate high-order correlations (Jovanović and Rotter 2016;Ocker et al. 2017). It would be interesting to consider the highorder statistics of stochastic bursting in the framework of the present model.
From the viewpoint of applications, we would like to mention that correlation functions of the spike trains are an important indicator of neural code in computational neuroscience (De La Rocha et al. 2007;Doiron et al. 2016;Nirenberg and Latham 2003). Our analytic results could be potentially compared with observed correlations and spectra, shedding light on the origin of observed correlations. In this context, it appears especially helpful that the basic expressions describing the spectral properties like (12), (14) and (25) contain only a few parameters. Therefore, experimentally observed spectra could be easily fitted to the derived analytic forms.