Controlling noise-induced behavior of excitable networks

The paper demonstrates the possibility to control the collective behavior of a large network of excitable stochastic units, in which oscillations are induced merely by external random input. Each network element is represented by the FitzHugh–Nagumo system under the influence of noise, and the elements are coupled through the mean field. As known previously, the collective behavior of units in such a network can range from synchronous to non-synchronous spiking with a variety of states in between. We apply the Pyragas delayed feedback to the mean field of the network and demonstrate that this technique is capable of suppressing or weakening the collective synchrony, or of inducing the synchrony where it was absent. On the plane of control parameters we indicate the areas where suppression of synchrony is achieved. To explain the numerical observations on a qualitative level, we use the semi-analytic approach based on the cumulant expansion of the distribution density within Gaussian approximation. We perform bifurcation analysis of the obtained cumulant equations with delay and demonstrate that the regions of stability of its steady state have qualitatively the same structure as the regions of synchrony suppression of the original stochastic equations. We also demonstrate the delay-induced multistability in the stochastic network. These results are relevant to the control of unwanted behavior in neural networks.


Introduction
An excitable system is a system with input and output, and its characteristic feature is that it demonstrates two drastically different kinds of response to different kinds of input: when input is lower than a certain threshold, the excitable unit demonstrates only small oscillations around the steady state, but if the input exceeds the threshold, the unit responds with a high spike whose duration and shape are almost independent of the shape or duration of the input signal. A single excitable unit is widely used as a simplified model of a typical neuron. Neural networks are responsible for information processing in animals (including humans) [1], and similar structures exist even in plants [2]. Because of that, their study and modeling attracts a lot of research attention. Networks of excitable units are popular models of the biological and artificial neural networks.
Depending on the properties of a single unit entering the network, and on the kind of coupling between them, the neural networks can display a rich variety of behaviors, but common to most of them are pattern formation and synchronization phenomena. The collective behavior of biological neural networks and of their models can be characterized by parameters such as timescale, regularity and synchronization strength. It is believed that synchronization in biological neural networks enhances information processing [3,4], but it is also linked to neural disorders like epilepsy and Parkinson's disease [5]- [8]. Therefore, depending on the circumstances one might wish to strengthen synchronization or to destroy it, or to change other parameters of the behavior. The desired control method should apply an intelligent weak perturbation capable of adjusting the behavior of the network gently, but efficiently. With this, it should not require any change in the connections or in the units within the network, and ideally not even detailed knowledge of the system structure. Finally, a good control tool should allow application on a macroscopic scale rather than on a scale of individual units.
Recently, a number of methods have been proposed for the suppression of synchrony in the arrays of coupled oscillators in which oscillations are self-sustained [9]- [11]. In this paper, we consider the possibility to manipulate the behavior of a large stochastic network of coupled excitable units, each demonstrating noise-induced spiking and unable to oscillate without external noise.
It was earlier proposed to control oscillations induced merely by random noise in a single [12]- [14], or in two mutually coupled [15] excitable units by means of delayed feedback [16,17] that was originally developed to control deterministic chaotic dynamics. The controlling force is proportional to a difference between the current system state and its state τ time units before. The main advantage of this technique is that it does not add new, or change the positions of the old, steady states, meaning that the disruption to the system is minimal. Also, for deterministic chaos it can be non-invasive. Indeed, a skeleton of a chaotic attractor consists of a countable set of unstable periodic orbits (UPOs), and if τ is equal to the period of one of them, the orbit may become stable. Also, it was recently shown that the delayed feedback can suppress any oscillations in a deterministic system altogether [18,19]. In both the latter cases, the controlling force vanishes, which is of course advantageous particularly in biological applications. However, when oscillations are induced by noise, there are no embedded periodic orbits, and although the control force can be minimized at optimal τ and strength, it never vanishes completely [12]- [15].
In the present work, we probe the delayed feedback as a possible tool to control the collective behavior of a large stochastic network of excitable units coupled through the mean field. This is motivated by the fact that, if the size of such a network tends to infinity, although each element behaves randomly, the mean field and other ensemble-averaged characteristics behave deterministically. Moreover, it was recently shown that the mean field of such a network can demonstrate deterministic chaos of various forms [20]. In this case, application of the delayed feedback to the mean field might produce effects similar to those observed in deterministically chaotic systems, implying the possibility of the desirable non-invasive control.
In the present work, we first reveal the effects produced by the delayed feedback in the stochastic excitable network numerically. We show that the feedback is capable of suppressing or weakening the collective synchrony in the network, and outline the areas on the plane of the feedback parameters where this effect takes place. Then we make an attempt to explain the numerically obtained results on a qualitative level. For this purpose, we employ the cumulant expansion of the probability density of the network, using Gaussian approximation, and derive the cumulant equations of the system with the feedback. We compare the results of analysis of the stochastic and cumulant equations and give evidence that the cumulant equations can explain the main features of the control, although only qualitatively.
Section 2 describes the model used and its collective behavior in the absence of delayed feedback. Section 3 demonstrates the numerically revealed effects of delayed feedback control applied to the network. Section 4 considers the dynamics of cumulant equations with delayed feedback. In section 5, the results of stochastic simulation and of the analysis of cumulant equations are compared and discussed.

Collective motion in the stochastic excitable network
We consider a network of excitable units coupled through the mean field, when each element is coupled to all other elements in the network with the same strength. This is probably the simplest form of coupling and represents a rough approximation of coupling in a real neural network. Each element number i of the network is modeled as a FitzHugh-Nagumo system influenced by noise, which is a paradigmatic model of a single excitable unit [21] namely, Here, ξ i (t) is Gaussian white noise with zero mean and unity variance, and the noise sources in different elements are uncorrelated, i.e. ξ i (t)ξ j (t + s) = δ i j δ(s); T is the noise intensity that is taken to be the same in all units. M X (t) is the averaged input from the rest of the network (mean field) and γ is the strength of coupling to the mean field. To ensure that uncoupled units cannot demonstrate repetitive spiking without external random perturbation, the value of |a| has to be larger than 1. Also, to enable separation of the xand y-timescales, ε has to be a small positive number, 0 < ε 1. In order to allow straightforward comparison with earlier works [20,22] where the networks of this type were considered (without control), we choose ε = 0.01 and a = 1.05. Note that a single spiking element placed inside the network will continue spiking, but the manner of its spiking may change considerably. As shown in [20]- [23], equations (1) can demonstrate various types of collective behavior depending on the level of noise and the strength of global coupling. Namely, different neurons can spike independently of each other, or spike mostly together, or be in any state in between these two extremes, i.e. spike simultaneously only occasionally. With this, the oscillations of the mean field can be of two main types: small subthreshold oscillations (no synchronization) and superthreshold spiking with large amplitude (synchronization). In their turn, both subthreshold and superthreshold oscillations can be either periodic, or irregular. As shown in [20], periodic spiking of the mean field corresponds to the largest degree of synchrony in the network, while chaotic spiking reflects weaker synchronization when the units spike together only at certain times. Also, the frequency of the mean field spiking is the largest when this spiking is periodic, i.e. when the network is maximally synchronized, as compared to the chaotic spiking.
The mean field obtained by numerical simulation of equations (1) is illustrated in figure 1. The top panel shows variance σ 2 X of M X (t) as a function of coupling strength γ at a fixed value of noise intensity T = 3.1 × 10 −4 for N = 100 (solid line) and for N = 5000 (circles). One can see that in both cases there is a range of values of γ in which synchronization occurs. When the number of coupled units is relatively small (N = 100), the transition from synchronous to non-synchronous behavior is gradual and is accompanied by the smooth decrease of σ 2 X to zero. However, in a considerably larger network (N = 5000) this transition is more abrupt. The behavior of networks of different sizes with otherwise the same parameters is compared in figure 1 (lower panels) where phase portraits and realizations of the mean field are given.
Two different kinds of irregular oscillations of the mean field for 10 000 oscillators are illustrated in figures 2(a) and (b) and 5(a) and (b). Namely, in figure 2(a) the mean field is spiking in a manner that reminds of chaotic oscillations, and this behavior will be further referred to as chaotic spiking (it was called intermittent spiking in [20]). In figure 5(a), the mean field is oscillating in a chaotic manner in the close vicinity of the fixed point. The former case corresponds to weak synchronization in the network, while the latter case to the absence of collective synchrony. In the next section these states will be taken as reference states, and we will find out if it is possible to change them by inducing, strengthening or weakening the existing synchrony.

Delayed feedback control of the collective behavior
In a real network it is natural to suppose that noise level and strength of coupling are likely to be determined by the environment and by the structure of the system, and it would not be easy to change their values without serious intervention with the network structure. However, the kind of behavior the network is displaying might not be satisfactory for some reason. For example, in humans the healthy state of the neural network of the brain is associated with non-synchronous oscillations of individual neurons. However, conditions like epilepsy, Parkinson's disease and essential tremor [5]- [8] are often associated with synchronization, which in this situation would have a detrimental effect to the functioning of the organism and needs elimination. On the other hand, synchronization is sometimes thought to enhance information processing by a group of neurons [24,25], and might therefore be desirable. In this context, the problem arises to control the properties of collective behavior of a neural network by some simple method requiring minimal invasion into the system. In the theory of control of deterministically chaotic oscillations a method is known which satisfies some of the requirements of the given task. It is the method of delayed feedback control originally proposed by Pyragas [16] with the aim of turning deterministically chaotic behavior into a periodic one. One registers an experimental signal s(t) coming from the chaotically oscillating system and forms a feedback force F(t) as follows: where s(t) is some experimental signal, K is the strength of the feedback and τ is time delay. It is well known that a chaotic attractor has a skeleton formed by a countable set of UPOs with different periods. If τ is equal to the period T of a certain UPO, then the given orbit may become stable and periodic oscillations with period T will be observed in the experiment. Remarkably, after all transients die out, the force F(t) vanishes and this control appears non-invasive, which is particularly appealing when it comes to medical and biological applications. However, there is another fact about the delayed feedback control which is of particular interest to us in the context of the given problem. In [19], the effect of the control force (2), on a deterministically chaotic system was studied for a large range of parameters K and τ . It was shown that on the (K , τ ) plane there are areas where all oscillations stop, and the feedback force F(t) tends to zero. These results are consistent with the earlier reported results on the feedback control of a periodic system [18].
The meaning of these results for a large stochastic network is as follows. Although each element of the network behaves stochastically, the probability density distributions (PDDs), the mean field and other ensemble averaged characteristics of the network are deterministic values, typically changing in time. As mentioned in section 2, the mean field can demonstrate a wide range of behaviors, including no oscillations, small periodic and chaotic oscillations, and large chaotic or periodic spiking-all of these regimes behaving deterministically. We hypothesize that if such a network is subjected to the delayed feedback control, in which the reference variable is the mean field, the ensemble averaged variables of the system should respond like the variables of a deterministic system. In particular, we can expect the suppression of oscillations of the mean field for some values of the feedback parameters, and this would imply desynchronization of the network.
Some other features of the delayed feedback control of a deterministic system include an essentially multi-leaf structure of the bifurcation diagram on the plane (K , τ ) [18,19], meaning that at the same values K , τ different regimes can be observed depending on the initial conditions. We can expect the same features to appear in the controlled stochastic network as well.
We introduce the delayed feedback in the equations for y-variable by analogy with [12]- [14] where the control of a single stochastic FitzHugh-Nagumo system was considered: Here, M Y is the mean field at the current time moment and M Y (t − τ ) is the delayed version of it. The feedback effects illustrated in this paper were found to be qualitatively the same when the feedback was introduced through M X in the first equation: the main difference was the feedback strength K required to produce similar changes. We first consider the effect of delayed feedback on the network equations (3) in the initial state of partial synchronization when the mean field demonstrated chaotic spiking (T = 0.000 28, γ = 0.1, N = 10 000, figures 2(a)). Here, and in other similar figures the length of each simulated realization was 60 000 time units. The map of regimes on the plane of control parameters τ and K is shown in figure 3, in which the mean spiking frequency is shown by color code. Black areas denote zero spiking frequency, which means the absence of spiking during the observation time and, consequently, the absence of synchronization in the network. One can see that synchronization can be suppressed even with a very small value of K if τ is chosen appropriately, which makes the control almost non-intrusive even before the transients die out. However, areas above the upper borderline of black areas correspond to the higher spiking frequency than in the initial state and therefore to stronger synchronization in the network, as shown in [20].
This diagram also illustrates multistability induced by delayed feedback in the network (see lower parts of the figure). Namely, both lower figures were obtained as τ was changed with a small step, and the final state of the system from the previous simulation was used as the initial state for the next calculation, which is an analogue of continuation procedure. The left figure was obtained as τ was increasing, and the right figure as τ was decreasing. One can see that these figures are different, in that the upper-right border of the suppression area has two different locations. Multistability occurs between the left-hand line C of some bifurcation (possibly non-local), as a result of which the current spiking regime vanishes, and the righthand line HB of Andronov-Hopf bifurcation of the steady state. Both lines are shown in the lower panels of figure 3 by continuous gray (green online) curves. This means that the system demonstrates bistability: at the same values of K and τ , different initial conditions can lead to either synchronous, or to non-synchronous states.
The map of regimes in figure 3 is in an excellent agreement with the earlier revealed bifurcation diagrams of deterministic systems with the feedback [18,19], at least in the areas where the steady state of the mean field is stabilized.
Two points on the map of regimes are illustrated in figures 2(c)-(f). Namely, in figures 2(c) and (d) K = 1, τ = 0.8 (which is outside the range of the diagram but in the region of enhanced  figure 4, where a few characteristics of the mean field oscillations are given as functions of τ . The convenient quantities that describe the essential properties of the network behavior are (i) the mean interspike interval (ISI) T of the spiking mean field M X , that defines the timescale of the synchronized oscillations; the standard deviation σ T of the ISI that characterizes the irregularity of spiking (the larger the σ T , the more irregular the spiking is); mean amplitude A of spiking that allows one to distinguish between the regimes with and without spiking; and the standard deviation σ X of the mean field M X that allows one to observe the transitions between the oscillations of different amplitudes. All these characteristics are given in figure 4, one can notice that variation of τ generally leads to the increase of both the period of oscillations (a) and the standard deviation of ISIs (b), which means less frequent and less regular spiking of the mean field. This, in its turn, implies weakening of synchronization by the feedback. Within certain ranges of τ (approximately from 0.25 to 1.2) the period is infinitely large, which, together with small σ X (d) is an indication of the complete loss of synchronization. Where synchronization is present, the mean field is spiking with approximately the same amplitude (c). Importantly, the increase of the mean period appears much more drastic than in the case of a single excitable unit controlled by the same method (compare with [12,13]). Next we consider the effect of delayed feedback on the initially non-synchronous network equations (3), when the mean field demonstrated chaotic subthreshold oscillations (T = 0.000 27, γ = 0.1, N = 10 000, figures 5(a)). The respective map of regimes on the plane (K , τ ) is shown in figure 6. One can see that the feedback can induce high-frequency spiking of the mean field, and thus strong synchronization in the network, at sufficiently large values of K . Also, there is an area of weak synchronization at smaller K , which is schematically indicated as a shaded region (green online) that is bounded by the black region of no spiking from below, and the region of high-frequency spiking from above. Two points on the map of regimes are illustrated in figures 5(c)-(f). Namely, figures 5(c) and (d) illustrate the weakly synchronous regime at K = 0.4, τ = 4. Figures 5(e) and (f) illustrate the strongly synchronous regime at K = 1, τ = 0.8 (this point is outside the diagram in figure 6), when the mean field spikes periodically.
The cut-off of the diagram along K = 0.1 is illustrated in figure 7 that shows the same characteristics as figure 4. As seen from the figure, the feedback is capable of inducing synchronization at certain values of τ (approximately from 1.2 to 2, see gray area in figures 7(a)-(d)). The average period T of the mean field spiking (a), as well as the standard deviation σ T of the ISI (b), has a minimum in the center of synchronization region which corresponds to the largest degree of synchrony. Inside synchronization domain, the spiking amplitude A remains constant within numerical accuracy. At the same time the standard deviation σ X of the mean field demonstrates a smooth maximum in the center of synchronization domain which indicates the highest spiking frequency and thus strongest synchronization. Finally, we examine two more initial states of the network at such values of γ = 0.1 and T when it is strongly synchronized. At T = 0.0003 the network is in a strongly synchronized state near the borderline of synchronization region [20], and the mean field is spiking periodically. Application of delayed feedback results in the map of regimes shown in figure 8(a). One can see that the area, where the synchrony is suppressed, is shrinked as compared to the case when initially the system was more weakly synchronized (compare with figure 6). When the system is deeper inside the synchronization region at T = 0.0005, delayed feedback cannot suppress synchronization completely (figure 8(a)), but can weaken synchronization, as indicated by the lower values of spiking frequency at nonzero values of K and τ .
The qualitative explanation of the observed effects of the delayed feedback will be attempted in the next section in terms of the dynamics of the cumulant equations.

Cumulant equations with delayed feedback
In the previous section, we reported the effects induced by delayed feedback in the stochastic network of excitable units modeled by equations (3). In particular, it was found that in spite of the random nature of oscillations in the network, the effects induced by the feedback are very similar to those in the deterministic systems. In this section, we will try to provide an approximate deterministic description of the stochastic network and thus to explain the effects induced by the feedback by analogy with deterministic systems.
The collective state of N stochastic units at a time moment t can be rigorously described by the joint PDD p 2N (x 1 , y 1 , x 2 , y 2 , . . . , x N , y N , t) of all system variables, which is a deterministic function of time. However, this method is inconvenient since p 2N is a function of 2N + 1 variables, where N is usually a large number. In order to simplify the description, one can assume that oscillations in different units of the network are uncorrelated to a certain degree, and that there are infinitely many coupled units, i.e. N → ∞. This is consistent with the approximation of molecular chaos often used in statistical physics in order to make the calculations tractable. Then one can represent p 2N as a product of N identical two-dimensional (2D) PDDs p 2 (x, y, t). The description of the network will then be reduced to the description of p 2 , which is a deterministic function of only three variables. This approximation was successfully used for noisy self-sustained oscillators [26], bistable units [27]- [29], phase oscillators [22,30,31], and excitable FitzHugh-Nagumo units [20,22,23] coupled via the mean field. In [32] stochastic bifurcations in coupled stochastic FitzHugh-Nagumo units were studied when N was allowed to be large, but finite.
The time evolution of the collective state of the network can be described by a Fokker-Planck equation for p 2 . However, a partial differential equation with delay is generally quite difficult for the analysis. Alternatively one can derive a system of equations governing the time evolution of cumulants of p 2 . A discrete set of cumulant first-order differential equations might be more convenient for the analysis than partial differential equations, but an inconvenience arises from the fact that generally there will be infinitely many cumulants. At this stage, some additional assumptions about p 2 might help to reduce the number of nonzero cumulants considerably and to make the analysis tractable.
The time evolution of p 2 is illustrated in figure 9, where several snapshots are shown at different time moments corresponding to different stages of oscillations of the mean field at T = 0.000 28, γ = 0.1, K = 0.1 and τ = 1.5. The function p 2 is numerically estimated from N = 10 000 units. The mean field is indicated by a filled circle. It is obvious already from this figure that p 2 has quite a complicated shape. To appreciate the complexity and non-smoothness of its shape, two 1D densities p X 1 (x, t) are shown in figure 10 for an uncoupled network (a) and for a network coupled with γ = 0.1 without delayed feedback. One can see that in a coupled network p X 1 has two maxima most of the time, which are sharp and non-equal, and the whole density is very non-symmetric. We would like to use an approximation of the p 2 by cumulants up to a certain order. Namely, given that there exists a system of infinitely many cumulants describing the density p 2 exactly, we could truncate this system at some stage, assuming that the sufficiently highorder cumulants tend to zero. The plots for the cumulants up to the sixth order are shown in figure 11. One can see that even the sixth-order cumulants do not show the tendency to go to zero. With this, the number and complexity of cumulant equations increase considerably with the order of cumulant expansion. For example, if we use cumulants up to the fifth order, the number of equations is 20 and the number of terms in them is just under 300. Such a complex nonlinear high-dimensional system is particularly difficult for the analysis if it includes delay terms.
With this, if the use of higher-order cumulants is not expected to provide a considerably better accuracy of the description of the function as illustrated in figure 11, it does not seem reasonable to make the problem so difficult for the purposes of qualitative analysis. One can well use a simpler approximation. The simplest approximation would be Gaussian and it assumes that p 2 can be approximately described as a Gaussian function of two variables. Then p X 1 (x, t) and p Y 1 (y, t) are approximately Gaussian functions, too. Gaussian distribution can be characterized by only five nonzero cumulants: mean values m X and m Y of variables x and y, respectively, their variances D X and D Y , and their cross-variance D X Y . If the Gaussian approximation is valid, the collective dynamics of the stochastic network of excitable units equations (3) can be described by a closed system of five differential equations for cumulants.
In [20,22] a Gaussian approximation was used to construct the cumulant equations for system (1) using the method proposed in [27] and improved in [33,34]. The validity of this approximation was verified by comparing the dynamics of cumulant equations with the dynamics of the mean field obtained by numerically solving the stochastic differential equations (SDEs) (1) at large N . It was found that the Gaussian approximation was capable of capturing the most essential features of the behavior of the stochastic network, at least on a qualitative level. In particular, the cumulant equations could demonstrate chaotic or periodic, sub-or superthreshold oscillations, just like the mean field in the original equations (1). Also, on the plane of parameters (T, γ ) similar transitions were observed for the mean field estimated from SDEs and from the approximate cumulant equations. These were the period-doubling bifurcations, saddle-node bifurcations, and some other, possibly non-local, bifurcations of the steady states and periodic solutions.
In this work we follow the example of [20,22] and exploit a Gaussian approximation to write down the cumulant equations for a network of FitzHugh-Nagumo oscillators equations (3) with delayed feedback. The appropriate averaging of stochastic equations with delay (3) was made using the techniques of [27,33,34] and the system of the following cumulant equations was obtained: We considered equations (4) at two sets of parameters (T, γ ), at which the system without feedback (K = 0) was in two distinct chaotic states: chaotic spiking (figures 12(a) and (b)) and subthreshold chaotic oscillations (figures 13(a) and (b)). The first state corresponds to the weakly synchronous network, and the second state to the non-synchronous one. For each regime above, we performed bifurcation analysis of equations (4) in the plane (τ, K ) with the help of a continuation technique using the free software DDE-BIFTOOL [35]. The resulting bifurcation diagrams are shown in figure 14. On the solid lines Andronov-Hopf bifurcation of the steady state takes place, and in gray areas the steady state is stable. In the white areas the oscillations of the mean field are subthreshold, either periodic or chaotic. The green points show the areas in which the mean field performs some kind of spiking, either periodic or chaotic. These areas were found by numerical simulation of cumulant equations rather than by continuation.
First, compare figure 14(a) with figure 3. The black areas of figure 3 correspond to white and gray areas of figure 14(a). We can see that for small K these areas have similar shapes, although at large K the similarity vanishes: in the stochastic equations the suppression regions are bounded from above (at K around 6 or less), while in the cumulant equations these regions show no tendency for becoming bounded from above in the range of K considered. Another  figure 14(a). From the earlier works on the analysis of deterministic systems with delayed feedback [18,19] it follows that the centers of the suppression regions are determined by the eigenperiods (2π divided by the imaginary part of the eigenvalue associated with Andronov-Hopf bifurcation) of the steady states without feedback. For example, the first area is formed around τ close to half of the eigenperiodτ . The same regularity is observed in the cumulant equations with the eigenperiod of the steady state being τ ≈ 0.72. Comparing to the SDEs, τ in the bifurcation diagram appears to be rescaled by a factor of ≈ 2, which can be explained by the inaccuracy of the Gaussian approximation, and also by the fact that the values of T in the stochastic and cumulant equations were taken to be different in order to reproduce qualitatively the same behavior. Next, compare the figure 14(b) with figure 6. At small and moderate values of K the shapes of the black regions in the latter match qualitatively the shapes of the white and black regions of the former. Also, at τ ≈ 0.65 there exists a region of spiking, which is comparable with the region of mean field spiking in figure 6 around τ ≈ 1.5. Again, at larger values of K the two diagrams do not match the ones for SDEs.
As predicted by [19], introduction of delayed feedback can induce multistability in the system which is illustrated in figures 15 and 16. This means that depending on the initial conditions, at exactly the same values of system parameters the system can demonstrate one of two (or of more) regimes. Interestingly, in both cases multistability occurs between two different spiking regimes: periodic and chaotic. The phase portraits of the two regimes do not look very different from each other, but the power spectra indicate the distinction between them clearly: they are discrete (up to numerical accuracy) in the case of periodic spiking and continuous with additional frequencies in the case of chaotic oscillations. The cumulant equations even of the second order, that approximately describe the large stochastic excitable network with delayed feedback, appear capable of capturing the most essential features of the map of regimes on the plane of the feedback parameters K and τ . Namely, they reproduce qualitatively the shape of the areas in which there is no synchronization, and therefore no mean-field spiking, at least for small and moderate values of K . There is no quantitative agreement, however, which is explained by the expected inaccuracy of the Gaussian approximation.

Summary
In this paper, we considered the possibility to employ delayed feedback in order to control the collective behavior of a large network of globally coupled stochastic excitable units, each modeling an excitable neuron in which spiking is induced merely by external random perturbation. It was essential that the delayed feedback was applied through the mean field. The action of control was initially studied by means of numerical simulation of the large system of stochastic differential equations with delay. A large range of feedback parameters was scanned, and the maps of regimes were obtained for a few distinct initial states of the network. Namely, it was shown that where the network was originally weakly synchronized, the feedback was capable of destroying synchronization even with a small value of K . Where the network was initially strongly synchronous, the feedback was shown to be capable of weakening synchronization.
On the other hand, where there was no synchronization from the beginning, the feedback could induce the latter at certain parameter values. In both cases, it was possible to adjust the network macro-parameters like timescales and the amount of coherence in the mean field. The explanation of the observed phenomena on a qualitative level was done by formulating and analyzing the simplified cumulant equations of the second order, which were reasonably tractable and allowed for bifurcational analysis using the tools developed for deterministic systems. In particular, it was shown that chaotic cumulant equations obey the same regularities that were previously found in other deterministically chaotic systems subjected to delayed feedback control. The main practical feature of the system being controlled is that there is quite a large area in the plane of the feedback parameters within which a fixed point is stabilized. In terms of the original network, this means desynchronization of coupled oscillations.
The results of this study are relevant to the task of non-intrusive control of populations of neurons, and in particular to the elimination of the unwanted synchronization that is thought to be linked to neural disorders like epilepsy and Parkinson's disease.