{'fraction_non_alphanumeric': 0.03811216429699842, 'fraction_numerical': 0.009281200631911532, 'mean_word_length': 4.686002994011976, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamically-Coupled Oscillators -- Cooperative Behavior via Dynamical Interaction --

We propose a theoretical framework to study the cooperative behavior of dynamically coupled oscillators (DCOs) that possess dynamical interactions. Then, to understand synchronization phenomena in networks of interneurons which possess inhibitory interactions, we propose a DCO model with dynamics of interactions that tend to cause 180-degree phase lags. Employing an approach developed here, we demonstrate that although our model displays synchronization at high frequencies, it does not exhibit synchronization at low frequencies because this dynamical interaction does not cause a phase lag sufficiently large to cancel the effect of the inhibition. We interpret the disappearance of synchronization in our model with decreasing frequency as describing the breakdown of synchronization in the interneuron network of the CA1 area below the critical frequency of 20 Hz.

Many studies of cooperative behavior in many-body systems, of which one prototype is represented by spin systems, have been based on the assumption that the dynamics of interactions can be ignored. However, in the modeling of many biological systems, the inclusion of interaction dynamics is necessary to obtain a faithful description. In much of the research that has been done, systems with such slow interaction have been mimicked by a so-called multi-compartment model which corresponds to a discrete version of a distributed parameter system. By using such a multi-compartment model, we can accurately estimate the dynamical behavior of a single unit. However, it is impossible to carry out a large-scale numerical simulation with such models to understand the behavior of a large-scale system composed of tens of thousands of units. Even if it is possible, we cannot understand the essence of such phenomena. To address this problem, we need to reduce the degrees of freedom needed in a descriptive model.
Phase reduction is a powerful tool that helps us understand the oscillatory behavior of large-scale systems composed of such multi-compartment models. In previous studies of systems consisting of such multi-compartment models, phase descriptions of the entire systems including the dynamical interactions themselves were obtained numerically [1,2,3,4]. Thus, the effects of the interaction dynamics were numerically incorporated into only a phase-coupling function to establish the phase equation. Therefore, in general, the effects of the interaction dynamics alone cannot be extracted in such methods.
In this paper, we analyze systems of dynamically coupled oscillators (DCOs), and demonstrate how the dynamics of the interactions influence the cooperative behavior exhibited by the systems. Using the multi-scale perturbation method (MSPM) [5,6] to carry out a phase reduction, we are able to represent the nonlinear dynamical interaction in terms of only the amplitude and phase of the fundamental frequency component in the frequency response of the interaction. In control engineering, the amplitude and phase of the fundamental frequency component of the frequency response are referred to as the describing function.
The describing function used in the treatment of non-linear systems is an extension of the transfer function used in the treatment of linear systems. In contrast to previous work [1,2,3,4], the effects of the interaction dynamics are expressed by the describing function in our approach, and so it is easy to isolate and extract the effects of the interaction dynamics from the phase equation. Therefore, our approach can clearly elucidate the effect of the interaction dynamics on the cooperative behavior of entire systems. Our results imply the possibility of a new modeling approach for biological systems by which the describing function of the interaction is identified. The description proposed here corresponds more to higher order approximation than to conventional linear description (e.g., the auto-regressive (AR) model).
In the first half of this paper, we show in detail the process of analytical phase reduction for a general model with dynamical interaction. We treat a large population of Stuart-Landau (SL) oscillators coupled through nonlinear dynamical interaction. The SL oscillators have the essential structure of the Hopf-bifurcation, because its evolution equation can be derived from any non-linear oscillator system with the Hopf-bifurcation through perturbation expansion [6]. In the latter half of the paper, we apply our theoretical framework to the analysis of a network of interneurons. The analysis of an interneuron network is a good example to verify the applicability of our approach for understanding the neural cooperative behavior via dynamical interaction.
First, we consider a large population of SL oscillators weakly coupled through interactions which themselves are nonlinear dynamical systems. We analyze the following coupled system: Here w j is the state variable (a complex number) of the jth oscillator (with a total of N). Ω is the average natural frequency and ǫω j represents the deviation from the average for the jth oscillator randomly distributed with a density represented by g(ω). The terms multiplied by ǫ are considered perturbations, and thus, the quantity ǫ controls the magnitude of perturbations. When ǫ = 0, the system has an unstable fixed point at the origin and a stable limit-cycle orbit on the unit circle in the complex plane represented by where φ j is the phase of the jth oscillator. When ǫ = 0, the system is neutrally stable with respect to a perturbation in the form of a temporal shift while it conserves a fixed orbit.
Thus, φ j depends only on the initial condition. The quantity s jk in Eq. (1) represents the output from the dynamical system of the interaction expressed by where x jk represents the internal state of the interaction, F jk (·) is a nonlinear function determining its dynamics, and V jk (·) is an output function.
We now discuss the dynamical interaction described by Eq.
, the higher harmonics resulting from the nonlinearity of this equation are superimposed on s jk . Here we restrict our consideration to the case that the output s jk is a periodic function possessing the same period as the input w k . In this case, s jk can be expanded into the following Fourier form:  (1) can be represented as where we write φ j as a function of ǫt to make explicit its slow evolution in time, and ǫu j is the first order change in w j introduced by the perturbation. If we substitute Eq. (4) into Eqs. (1) and (3) and expand the resulting expression about ǫ = 0, the O(ǫ) equation Here L φ is the linear operator corresponding to Eq. (1) linearized about the periodic solution (2) in the case ǫ = 0. It is given by where the overline denotes complex conjugation. All of the eigenvalues of L φ are non-positive, since the solution Φ(t) is stable. Note that there exists the eigenfunction Φ ′ (t) of L φ with eigenvalue 0. This eigenfunction corresponds to an infinitesimal temporal shift because We assume there exists no other eigenfunction with eigenvalue 0 in the space of the periodic functions, so that kerL φ = span{Φ ′ (t)} = span{iΦ(t)}. This assumption is equivalent to that of the orbital stability of Φ(t).
We define the inner product of two 2π Ω -periodic complex functions u(t) and v(t) as 2π Ω 0 dt (Re{u(t)}Re{v(t)} + Im{u(t)}Im{v(t)}). Then, we obtain the operator adjoint to L φ as From Fredholm's alternative [5], there is a so-called response function Φ * (t) that spans the kernel of L * φ in the space of periodic functions. Therefore, kerL * φ = span{Φ * (t)}. In this case, we explicitly obtain Φ * (t) = iΦ(t). Taking the inner product of iΦ(t) and Eq. (5), we where the quantity ǫt is treated as a constant over a single period, since φ j (ǫt) is driven by a weak perturbation. Here, the higher harmonics resulting from the nonlinearity of Eq. (3) cancel, because iΦ j , J (n) jk Φ n k = 0 (n > 1). Thus, we derive the phase equation describing the slow phase dynamics of the system, with A jk = |J where ω j represents a natural frequency of unit j at the level of the phase equation. As Eqs. (9) and (10)  When |m| = 0, the system is in a state of phase synchronization. In the thermodynamic limit, N → ∞, we obtain the following self-consistent equation relating the order parameters |m| andΩ: Here, to derive the order parameter equation, we have assumed the existence of one large cluster of oscillators phase-locked at frequencyΩ in the slow dynamics described by (9). In the original system, described by (1), the frequency of the synchronous cluster is Ω + ǫΩ.
Next, we apply our theoretical framework to the analysis of a network of interneurons as an example. High frequency oscillations, known as γ oscillations, have been observed in many areas of the brain. The synchronization of inhibitory interneurons is known to be the origin of γ oscillations [8]. Inhibitory interactions, which correspond to anti-ferromagnetic interactions in spin systems, introduce competition; for this reason, in general, such interactions tend to prevent synchronization. Therefore, the synchronization phenomena in networks of interneurons are non-trivial [2,9,10]. A series of physiological in vitro experiments inspired our study. One in particular was the experimental discovery that the interneuronal network in the CA1 area of the hippocampus displays synchronization at frequencies greater than or close to 20 Hz, but not at frequencies significantly below this value [11,12]. The discovery of this phenomenon provides evidence suggesting that cooperative behavior in networks of interneurons arises due to the dynamical nature of the interactions.
To understand the mechanism of the synchronization, we propose a DCO model with interaction dynamics that tend to cause 180-degree phase lags. The estimated connectivity between interneurons in the CA1 area is 10% [10]. The interneuron network in the CA1 area has dense connectivity, because the number of connections between interneurons is O (N). Therefore, we can capture properties of the interneuron network by analyzing the global coupling systems. As shown in Fig. 1(b), the dynamical interactions of neurons are modeled by a rectangular-pulse generator consisting of a threshold element and a secondorder lag system (SOLS): The parameters τ (1) jk and τ (2) jk are the time constants of the two cascade units of the first-order lag system that form the SOLS. In Eq. (12), the parameters G and H denote, respectively, the gain and the threshold of the threshold element, and are defined as G = 1/θ and H = cos θ (θ > 0). If the value of θ is small, the rectangular pulse generated by the threshold element is sharp. Note that Eq. (12) represents mutual interactions only between the real parts of the SL oscillators.
In the context of control engineering, the SOLS is regarded as one of the minimal models for realizing a 180-degree phase lag. As shown in Fig. 1(a), the electrical properties of somata and proximal dendrites can be approximated as the CR circuit. Near resting potential, we can obtain a linear approximation of their behavior as being that of a first-order lag system.
The response of the synaptic conductance is also displayed in Fig. 1(a). The evolution of the synaptic conductance is roughly reproduced by that of the first-order lag system. Thus, as shown in Fig. 1(a), the SOLS given by Eq. (12) can be considered as a cascade system consisting of these components. In addition, to model the generation of action potentials in a neuron when its membrane potential reaches the threshold, and furthermore, to demonstrate the applicability of our approach in analyses of non-linear systems, we have also introduced a threshold element into the dynamical interaction given in Eq. (12).
Here, the describing function of the threshold element is 2 sin θ/(πθ), which contributes only to the gain of Eq. (13). Figure 1(c) displays the Bode diagrams of the SOLS. As Fig. 1(c) reveals, the gain of the SOLS, which is a type of low-pass filter, is a decreasing function of frequency at high frequencies, while the phase of the SOLS converges to -180 degrees as the frequency increases. If the phase lag is 180 degrees, J jk is a positive real number, and the system thus effectively becomes a ferromagnetic oscillator system. Note that Eq. (12) represents the mutual interactions only between the real parts of the SL oscillators. Thus, A jk takes half the value it had in Eq. (10), and here A jk = 1 2 |J jk . In the numerical calculation we will now discuss, we chose the distribution of natural frequencies as g(ω) = (2πσ 2 ) −1/2 exp (−ω 2 /(2σ 2 )), where σ is the deviation of the natural frequencies. We begin with systems where the interactions are all equivalent, i.e., Ωτ jk = Ωτ , and hence A jk = A and ψ jk = ψ. In this case, Eq. (9) is equivalent to the Sakaguchi-Kuramoto model [7]. Figure 2(a) displays |m| as a function of σ in the case of the dynamical interaction given by Eq. (12), where the solid curves were obtained from the order parameter equation (11), and the data points with error bars represent results obtained using the fourth-order Runge-Kutta method with Eqs. (1) and (12). Next, we estimate the critical value of σ, representing the boundary between |m| = 0 and |m| = 0.  σ c is proportional to A [7]. Plotted in Fig. 2(b) is the critical deviation normalized to have a maximum value of 1. The widths of the synchronous regions decrease as the frequency increases because, as shown in Fig. 1(c), the SOLS is a low-pass filter whose gain rapidly decreases at higher frequencies. However, if G in Eq. (12) is sufficiently large, the gain of the interaction given by Eq. (12) can be made sufficiently large by the pulse generator to realize synchronization at high frequencies. However, as Fig. 2(b) reveals, our model does not exhibit synchronization at low frequencies because, as shown in Fig. 1(c), this dynamical interaction does not cause a phase lag sufficiently large to cancel the effect of the inhibition. We interpret the disappearance of synchronization in our model with decreasing frequency as describing the breakdown of synchronization in the interneuron network of the CA1 area below the critical frequency of 20 Hz [11,12]. Based on this correspondence, we conjecture that in an interneuron network, a phase lag resulting from the dynamical nature of the interactions might cancel the effect of inhibition at high frequencies, and through this cancellation, the system effectively becomes a ferromagnetic oscillator system.
By means of numerical simulations, as shown in Fig. 3, we have demonstrated through a scaling plot that the variation of a local field is O(1/ √ N ) in systems where a quenched random time constant is used for the interaction dynamics in Eq. (12). In the thermodynamical limit N → ∞, the system asymptotically approaches a system with global homogeneous interaction, even with a quenched random time constant of the interaction dynamics. Therefore, as the system size N increases, the system tends to synchronize through the cancellation