{'fraction_non_alphanumeric': 0.04083906915765323, 'fraction_numerical': 0.018551294657489347, 'mean_word_length': 4.321241716079526, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of delay-coupled excitable neural systems

We study the nonlinear dynamics of two delay-coupled neural systems each modelled by excitable dynamics of FitzHugh-Nagumo type and demonstrate that bistability between the stable fixed point and limit cycle oscillations occurs for sufficiently large delay times and coupling strength. As the mechanism for these delay-induced oscillations we identify a saddle-node bifurcation of limit cycles.


INTRODUCTION
The brain may be conceived as a dynamic network of coupled neurons [1,2,3]. These neurons are excitable units which can emit spikes or bursts of electrical signals. In order to describe the complicated interaction between billions of neurons in large neural networks, the neurons are often lumped into highly connected sub-networks or synchronized sub-ensembles. Such neural populations are usually spatially localized and contain both excitatory and inhibitory neurons [4].
The simplest model to display features of neural interaction consists of two coupled neural systems. Starting from this simplest network motif, larger networks can be built, and their effects may be studied. For example, starting from two interconnected reticular thalamic neurons with oscillatory behavior, it was shown in [5] how more complex dynamics emerges in ring networks with nearest neighbors and fully reciprocal connectivity, or in networks organized in a two-dimensional array with proximal connectivity and "dense proximal" coupling in which every neuron connects to all other neurons within some radius. In another example [6], a neural population was itself modeled as a small sub-network of excitable elements, to study hierarchically clustered organization of excitable elements in a network of networks.
Most studies represent a population as an effective oscillatory element that is then coupled with other populations to study synchronization and network effects [7,8,9]. In these studies, the basis for the emergence of complex network dynamics is the oscillatory behavior of neural elements, in spite of the fact that individual neural systems are usually in a stable steady state exhibiting excitable dynamics when perturbed. A reason is that self-sustained oscillatory behavior in the individual element is required to define synchronization of subsystems [10]. One way to resolve this is to add weak Gaussian white noise to each individual element to generate sparse Poisson-like irregular spiking patterns, as seen in real neurons [6,11,12].
In this report, we introduce a time-delayed coupling into the model of two excitable populations and demonstrate that the coupling-delay can induce sustained oscillations between the two subsystems. We find a regime of bistability between the stable fixed point and limit cycle oscillations for sufficiently large delay times τ and coupling strength C. As the mechanism for these delayinduced oscillations we identify a saddle-node bifurcation of limit cycles.
This result suggests the use of a compound system of two time-delayed coupled excitable elements as a minimal network motif to investigate oscillatory behavior in more complex networks.

MODEL
We examine the delayed linear symmetric coupling of two identical neural populations. Each population is represented by a simplified FitzHugh-Nagumo (FHN) system [13,14], which is widely used as a paradigmatic model of excitable systems [15].
The dynamical equations are given by: where the two subsystems (x 1 , y 1 ) and (x 2 , y 2 ) correspond to two neuron populations. The value of the excitability parameter a determines whether the subsystem is excitable (a > 1) or exhibits self-sustained periodic firing (a < 1). This is because at a = 1 the uncoupled systems exhibit a Hopf bifurcation, and the fixed point becomes an unstable focus for a < 1. ǫ is the timescale parameter that, if chosen to be much smaller than unity, results in fast activator variables x 1 , x 2 , and slow inhibitor variables y 1 , y 2 . Unless otherwise noted, we shall choose ǫ = 0.01 and a = 1.3 for numerical simulations and thus restrict our analysis to parameter values where each of the two subsystems exhibits excitability with a stable fixed point.
The interaction between the two neural populations is modelled as diffusive, i. e., the coupling vanishes if the variables x i are identical. The coupling strength C is taken to be symmetric for simplicity. Further, we assume a finite signal propagation speed between distant neural populations. This is incorporated using the delay time τ . Note that the coupling term has the form of a classical diffusion term, but has lost its diffusive character through the introduction of a propagation delay τ . More general delayed couplings are studied in [16].
Before we analyse Eq. (1), let us briefly give a biophysical interpretation of the hitherto abstract nature of the dynamical variables x i and y i . This is needed to relate our results to other work on population dynamics in neural networks, in particular to studies dealing with spiking rates and time delays [17,18,19,20]. These studies consider neural fields [21,22], that is, they extend over one or more spatial dimensions. While neural field models adopt the continuum limit of a network, we consider the minimal network motif, i. e., two discrete subsystems, but in both cases the biophysical reason to include delay is the same.
The generic biophysical interpretation of the FHN model is based on a single point-like neuron and was originally not derived from features of a neural population.
The variable x models fast changes of the electrical potential across the membrane (spikes), and y is related to the gating mechanism of membrane channels [13,14]. On the contrary, a neural population is generically described by a time coarse-grained mean-field model, in which the dynamical variables represent averaged spiking rates [4]. Notwithstanding, spiking rates in neural populations can exhibit both steady state values and relaxation oscillations. Both states are described by the FHN mechanism with a > 1 and a < 1, respectively, which justifies the FHN model for "point-like" populations. To avoid confusion about spikes vs. spiking rates, a model of uncoupled subsystems as defined by Eq. (1) with C = 0 is sometimes referred to as the Bonhoeffer-van der Pol model [23,24,25], for example in Ref. [26].

LINEAR STABILITY ANALYSIS
The unique fixed point of the system is symmetric and is given by Denoting for convenience x(t−τ ) ≡x, we can re-write system Eq. (1) as: This system can be linearized around the fixed point x * by setting x(t) = x * + δx(t): with the Jacobian matrices J * f andJ * g . The explicit form where u is an eigenvector of J * f implies This leads to the characteristic equation which can be factorized giving This transcendental equation has infinitely many complex solutions λ. Fig. 1 shows the real parts of λ for various values of C. As can be seen in Fig. 1 the real parts of all eigenvalues are negative throughout, i. e., the fixed point of the coupled system remains stable for all C. This can be shown analytically for a > 1 by demonstrating that no delay-induced Hopf bifurcation can occur. Substituting the ansatz λ = iω into Eq. (8) and separating into real and imaginary parts yields for the imaginary part This equation has no solution for a > 1 since |ξ| = a 2 − 1 + C > C, which proves that a Hopf bifurcation cannot occur.

DELAY-INDUCED OSCILLATIONS
The cooperative dynamics of delay-induced oscillations in excitable systems is inherently different from those of noise-induced oscillations. The introduction of noise terms induces sustained oscillations in each individual subsystem by continuously kicking these subsystems out of their respective rest states. Coupling then produces synchronisation effects between these individual oscillators [11,12]. If time-delayed feedback control [27] is applied locally to one of the subsystems, the stochastic synchronisation can be tuned by varying the timedelay. This is in line with other work where it was demonstrated that such time-delayed feedback can be used to control the coherence and the timescales of noiseinduced oscillations in a single FitzHugh-Nagumo system [28,29,30,31].
For delayed coupling the case is entirely different. Here the sustained oscillations are an effect of the cooperative dynamics. They are generated by the delayed interaction between two non-oscillating stable units, and are thus an emergent phenomenon of the compound system. The bifurcation parameters for delay-induced bifurcations are the coupling parameters C and τ .
For large coupling delay τ the oscillation is readily understood as the two units firing alternately, each spike initiated by the delayed signal of the remote system. Since the signal of one system is transmitted to the other and then back, the oscillation in one system must have a period T of approximately 2τ , and be phase shifted by T /2 with respect to the other system. This is visible in the time series (Fig. 2(a),(b)) and in the phase portraits of the activators (Fig. 2(c)) and inhibitors (Fig. 2 (d)).
Since the two subsystems have identical but T /2shifted profiles x i (t) and y i (t), the coupling terms C[x i (t−τ )−x j (t)] (i = j) would vanish for all t if the oscillation period T were precisely 2τ . Then the periodic orbit ( Fig. 3(b)) would also be a solution for the uncoupled subsystems (C = 0). This cannot be the case for the excitable regime of the FHN systems (a > 1) because periodic orbits do not exist there. Hence it follows that the oscillation period T must be of the form T = 2(τ + δ), where δ = 0 is the effective time shift responsible for the non-vanishing coupling term. The phase portraits of the activators x i in the planes (x 1 (t), x 2 (t − τ )) deviate from the bisector (Fig. 4(a)) due to the fact that δ = 0.
To obtain a clear picture of the timescales involved in the dynamics, we have computed the excursion times along the segments of the phase-space trajectory. The start and end points of the different segments of the trajectory (colored dots A to D) are defined to correspond to the time steps when the trajectory has left or entered a neighborhood ∆x 1 (here ∆x 1 = 0.01) of the stable branches of the x 1 -nullcline.
In the example shown in Fig. 3 (ǫ = 0.01, a = 1.3, C = 0.5, τ = 3) the value of the effective time shift is δ = 0.012. This value is about one third of the fast transition times between the stable branches of the x 1 -nullcline, ∆t 1 = 0.041 and ∆t 3 = 0.038, respectively. ∆t 1 is the rise-time of the spike, i.e., the time that elapses between leaving the fixed point (black dot A) and crossing the right stable branch of the x 1 -nullcline at B (green dot, Fig. 3(b)). ∆t 3 is the drop-time, i.e., the duration of the jump back from the right to the left stable branch of the x 1 -nullcline (between blue and magenta dots, C → D). The slow parts of the trajectory, occurring on the right (B → C) and left (D → A) stable branches of the x 1 -nullcline, have a duration ∆t 2 = 0.357 and ∆t 4 = 5.588, respectively. The total oscillation period is thus T = 6.024. In Fig. 5 the phase portraits of the first of the two delay-coupled subsystems are shown for different excitability parameters a and delay times τ . The top panel (a,b) corresponds to excitability parameters far from the Hopf bifurcation of the uncoupled system, which occurs at a = 1. The bottom panel (c,d) corresponds to values of a close to the Hopf bifurcation.
In the case of a = 1.3 and τ = 3 ( Fig. 5(a)), the oscillation period T = 6.024 is large enough for the two subsystems to nearly approach the fixed point before being perturbed again by the remote signal. Note that Fig. 5(a) corresponds to Fig. 3(b). If the delay time becomes much smaller, e. g., for τ = 0.8 ( Fig. 5 (b)), the excitatory spike of the other subsystem arrives while the first system is still in the refractory phase, so that it cannot complete the return D → A to the fixed point. The times spent in the different stages of the limit cycle for τ = 0.8 are ∆t 1 = 0.087, ∆t 2 = 0.111, ∆t 3 = 0.085, and ∆t 4 = 1.355, and the oscillation period is T = 1.637. The effective time shift δ = 0.018 for τ = 0.8 is larger than for τ = 3.
We now investigate the transition from large to small delay times τ close to the Hopf bifurcation at a = 1.05. Again we choose τ = 3 (Fig. 5(c)) and find T = 6.018 and δ = 0.009. The times spent in the different phases of the oscillation period T are ∆t 1 = 0.043, ∆t 2 = 0.329, ∆t 3 = 0.186, and ∆t 4 = 5.461. For the shorter delay time τ = 0.8 (Fig. 5(d)) we find T = 1.630, δ = 0.015, ∆t 1 = 0.072, ∆t 2 = 0.240, ∆t 3 = 0.070, and ∆t 4 = 1.249. We find the same pattern as far from the Hopf bifurcation. The effective time shift δ is smaller if τ becomes larger and the system can approach the fixed point. Furthermore, if the system is close to the Hopf bifurcation  and τ is sufficiently large so that the fixed point is very closely approached, δ becomes much smaller than far from the Hopf bifurcation, and tends to zero for a → 1.
The reason is that when the excitatory spike arrives in the first subsystem, there is a turn-on delay δ before the first subsystem emits a spike. This is because the trajectory has to cross the middle branch of the x 1 -nullcline due to this upstream impulse. This section of the trajectory becomes the smaller, the closer the fixed point A (at x 1 = −a) is to the minimum of the x 1 nullcline (at x 1 = −1). This also explains the origin of the time-shift δ.
Finally, we shall investigate the question whether the system exhibits bistability between the fixed point and the limit cycle oscillation for all values of τ . In Fig. 6 the regime of oscillations is shown in the parameter plane of the coupling strength C and coupling delay τ . The oscillation period is color coded. The boundary of this colored region is given by the minimum coupling delay τ min as a function of C. For large coupling strength τ min is almost independent of C; with decreasing C it sharply increases, and at some minimum C no oscillations exist at all. At the boundary, the oscillation sets in with finite frequency and amplitude as can be seen in the inset of Fig. 6 which shows a cut of the parameter plane at C = 0.8. The oscillation period increases linearly with τ . The mechanism that generates the oscillation is a saddle-node bifurcation of limit cycles (see Fig. 6 inset (b)), creating a pair of a stable and an unstable limit cycle. The unstable limit cycle can be visualized from transients starting at appropriate initial conditions (see   7). The trajectory initially remains near the unstable limit cycle, which separates the two attractor basins of the stable limit cycle and the stable fixed point, and then asymptotically approaches the stable fixed point.

CONCLUSION
We have shown that delayed coupling can induce periodic spiking in a compound system of two coupled neuron populations if the delay and the coupling strength are sufficiently large. Bistability of a fixed point and limit cycle oscillations occur even though the single excitable element displays only a stable fixed point. The two neural populations oscillate with a phase lag of π.