Reduction Theories Elucidate the Origins of Complex Biological Rhythms Generated by Interacting Delay-Induced Oscillations

Time delay is known to induce sustained oscillations in many biological systems such as electroencephalogram (EEG) activities and gene regulations. Furthermore, interactions among delay-induced oscillations can generate complex collective rhythms, which play important functional roles. However, due to their intrinsic infinite dimensionality, theoretical analysis of interacting delay-induced oscillations has been limited. Here, we show that the two primary methods for finite-dimensional limit cycles, namely, the center manifold reduction in the vicinity of the Hopf bifurcation and the phase reduction for weak interactions, can successfully be applied to interacting infinite-dimensional delay-induced oscillations. We systematically derive the complex Ginzburg-Landau equation and the phase equation without delay for general interaction networks. Based on the reduced low-dimensional equations, we demonstrate that diffusive (linearly attractive) coupling between a pair of delay-induced oscillations can exhibit nontrivial amplitude death and multimodal phase locking. Our analysis provides unique insights into experimentally observed EEG activities such as sudden transitions among different phase-locked states and occurrence of epileptic seizures.

Since we consider global coupling, we assume K jk = K and L jk = L for all j and k in Eq.
(2). The number of the oscillators used in the numerical simulations is N = 256. The mean frequency of the oscillators is chosen as the frequency obtained at the point A in Fig. 1, which we denote by Ω. The frequency of each oscillator is independently chosen around this Ω from a normal distribution with a fixed standard deviation σ. For all oscillators, we assume the same value of the control parameter μ, i.e., their distances from the bifurcation points are the same as that of the point A. The parameter sets (α j β j ) satisfying such conditions are then inversely determined from Eq. (3) and from the relation α = μ + A, while keeping other parameters fixed at the same values as those in the main text, i.e., γ = −2, δ = 0, = −10, and t 0 = 8.

A. Kuramoto transition in a population of delay-induced oscillations
When the mutual coupling is weak enough, oscillators with nonidentical frequencies behave incoherently and no macroscopic oscillations are observed. As the coupling strength is increased, the oscillators tend to synchronize with each other and macroscopic coherent rhythms emerge. The most well-known example is the macroscopic synchronization transition in systems of globally coupled oscillators investigated by Kuramoto [1]. Suppose a system of N coupled oscillators described by the following phase model: where φ j is the phase of the jth oscillator, ω j is the natural frequency, and K p is the coupling strength. The natural frequency ω j of each oscillator is independently drawn from a probability density g(ω), which is assumed to be symmetric and unimodal with a peak at ω = ω 0 . It can be shown that this system undergoes a macroscopic synchronization transition at K p = 2/(πg(ω 0 )) and exhibits collective rhythms.
We here consider the Kuramoto transition in a system of globally coupled delay-induced oscillations described by Eq. (1) in the main text. When the frequencies of the oscillators are narrowly distributed and the coupling strength is relatively weak, we can approximate Eq. (1) by the reduced phase equation (6) and treat the collective dynamics of the system analytically. We assume that the oscillators interact only through the x component, i.e., we assume K > 0 and L = 0, and vary K as the control parameter. Equation (6) can then be converted to the Kuramoto model given above. The critical coupling strength K c is predicted in Ref.
[1] to be To check this prediction, we performed numerical simulations of Eq. (1). The parameter σ was set at σ = 1.5 × 10 −3 . Figure S1(A) shows the parameter sets (α j , β j ) used in the simulations. Figure S1(B) shows the time series of some of the delay-induced oscillations for K < K c (top panel) and for K > K c (bottom panel). We can clearly see that the oscillations are mutually synchronized when K > K c , though their waveforms are not completely the same because their parameters are slightly different. Degree of the macroscopic synchronization is measured by the order parameter R = 1 N | N j=1 x j | , where ... denotes time averaging. If all oscillators are perfectly synchronized at the same phase,

B. Amplitude death in a population of delay-induced oscillations
When the frequency distribution of the oscillators is relatively wide, the amplitude death phenomenon, in which all oscillators are stabilized at their origin and stop oscillations, is expected. Matthews, Mirollo, and Strogatz [2] carefully investigated the collective dynamics of the CGL Eq. (5) and derived the conditions for the amplitude death. It can be considered a generalization of the amplitude death phenomenon for two-oscillator systems that we treated in the main text.
For simplicity, we consider a specific case that the ratio of the coupling strengths in the x andẋ components is kept constant, namely, we assume that K and L are given by with M being a control parameter. Equation (1) can then be reduced to the CGL equation (5) with a coupling term (Mμa/2N)(u k − u j ), which is the form investigated in Ref. [2].
It is predicted that, under a necessary condition the amplitude death occurs when the control parameter M exceeds 1.
We fix the parameter σ at σ = 1.5 × σ c for the numerical simulations, which results in the parameter sets (α j , β j ) shown in Fig. S1 namely, the amplitude death actually occurs in a population of delay-induced oscillations.