Abrupt Desynchronization and Extensive Multistability in Globally Coupled Oscillator Simplices

TL;DR

This paper investigates the complex collective dynamics of coupled phase oscillators with three-way interactions, revealing abrupt desynchronization transitions and extensive multistability, which are relevant for understanding systems like the human brain.

Contribution

It introduces a novel analysis of simplicial interactions in oscillator systems, uncovering unique phenomena such as continuum desynchronization and infinite multistability.

Findings

Discovered abrupt desynchronization transitions without synchronization counterparts.

Identified extensive multistability with infinitely many stable states.

Highlighted implications for systems like the human brain and social networks.

Abstract

Collective behavior in large ensembles of dynamical units with non-pairwise interactions may play an important role in several systems ranging from brain function to social networks. Despite recent work pointing to simplicial structure, i.e., higher-order interactions between three or more units at a time, their dynamical characteristics remain poorly understood. Here we present an analysis of the collective dynamics of such a simplicial system, namely coupled phase oscillators with three-way interactions. The simplicial structure gives rise to a number of novel phenomena, most notably a continuum of abrupt desynchronization transitions with no abrupt synchronization transition counterpart, as well as, extensive multistability whereby infinitely many stable partially synchronized states exist. Our analysis sheds light on the complexity that can arise in physical systems with simplicial interactions like the human brain and the role that simplicial interactions play in storing information.

Equations

$${\dot{\theta }}_i$$

$${\varphi }_i$$

$${\dot{\theta }}_i={\omega }_i+K{r}_1^2\sin [2({\psi }_1-{\theta }_i)].$$

$${\partial }_{t}a=-2i\omega a+K\left(z_1^{\ast 2}-z_1^2{a}^2\right).$$

$${\dot{z}}_2=2i{\omega }_0{z}_2-2\Delta z_2+K\left(z_1^2-z_1^{\ast 2}{z}_2^2\right).$$

$${\dot{r}}_2$$

$${\dot{\psi }}_2$$

$$f_{\text{locked}}(\theta ,\omega )=\eta \delta (\theta -{\theta }^{\ast }(\omega ))+(1-\eta )\delta (\theta -{\theta }^{\ast }(\omega )-\pi ),$$

$$f_{\text{drift}}(\theta ,\omega )=\frac{\sqrt{{\omega }^2-K^2{r}_1^4}}{2\pi \left[\omega +K{r}_1^2\sin (2{\psi }_1-2\theta )\right]}.$$

$$r_1=(2\eta -1){\int }_{-K{r}_1^2}^{K{r}_1^2}\sqrt{\frac{1+\sqrt{1-(\omega /K{r}_1^2{)}^2}}{2}}g(\omega )d\omega ,$$

$$r_2=\frac{-1+\sqrt{1+K^2{r}_1^4}}{K{r}_1^2}.$$

$$\sqrt{\frac{2r_2}{1-r_2^2}}=\sqrt{K}(2\eta -1)$$

$$\times {\int }_{-2r_2/(1-r_2^2)}^{2r_2/(1-r_2^2)}\sqrt{\frac{1+\sqrt{1-{\left[\omega (1-r_2^2)/2r_2\right]}^2}}{2}}g(\omega )d\omega .$$

$${\eta }_{\text{min}}(K)=\frac{\sqrt{K_c(1)}}{2\sqrt{K}}+\frac{1}{2}.$$

$$r_1^{\text{min}}(K)\approx \frac{1.2120}{\sqrt{K}}\text{and}{r}_2^{\text{min}}(K)\approx \mathrm{0.5290.}$$

$$K_c(\eta )\approx \frac{4.137}{(2\eta -1{)}^2}.$$

Full text

Research into the macroscopic dynamics of large ensembles of coupled oscillators have extended our understanding of natural and engineered systems ranging from cell cycles to power grids Strogatz2003 ; Pikovsky2003 ; Prindle2012Nature ; Rohden2012PRL ; Skardal2015SciAdv . However, with few exceptions (including Tanaka2011PRL ; Komarov2015PRE ; Bick2016Chaos ), little attention has been paid to the synchronization dynamics of coupled oscillator systems where interactions are not pair-wise, but rather $n$-way, with $n\geq 3$. Such interactions are called “simplicial”, where an $n$-simplex represents an interaction between $n+1$ units, so $2$-simplices describe three-way interactions, etc Salnikov2019EJP . Recent advances suggest that simplicial interactions may be vital in general oscillator systems Ashwin2016PhysD ; Gin2018 ; Gin2019 and may play an important role in brain dynamics Petri2014Interface ; Giusti2016JCN ; Sizemore2018JCN and other complex systems phenomena such as, the dynamics of collaborations Patania2017EPJDS or social contagion Iacopini2019 . In particular, interactions in 2-simplices (named holes or cavities) are important because they can describe correlations in neuronal spiking activity (that can be mapped to phase oscillators Politi2015 ) in the brain Reimann2017 providing a missing link between structure and function. In fact, coupled oscillator systems that display clustering and multi-branch entrainment have been shown to be useful models for memory and information storage Daido1996PRL ; Hoppensteadt1999PRL ; Fell2011Nature ; Hipp2012Nature ; Ashwin2007SIAM ; Skardal2011PRE ; Komarov2013PRL . Despite these findings, the general collective dynamics of coupled oscillator simplices and their utility in storing information are poorly understood.

In this work we study large coupled oscillator simplicial complexes, considering the impact of $2$-simplices, i.e., three-way interactions, on collective behavior. Specifically, we consider the 2- and 1-simplex multilayer system given by

[TABLE]

The dynamics in the $\theta$-layer are the natural generalization of the classical Kuramoto model Kuramoto1984 with $2$-simplex interactions (namely, coupling is sinusoidal and diffusive), where $\theta_{i}$ represents the phase of oscillator $i$ with $i=1,\dots,N$, $\omega_{i}$ is its natural frequency which is assumed to be drawn from the distribution $g(\omega)$, and $K$ is the global coupling strength. The $\theta$-oscillators then drive a population of $\phi$-oscillators in another layer, for which a one-to-one correspondence exists following a multiplex structure DeDomenico2013PRX . Specifically, the $\phi$-oscillators evolve subject to the classical Kuramoto model (with natural frequencies $\nu_{i}$ and coupling strength $\kappa$), but with an additional driving term with strength $d$.

While numerical investigations of systems with non-pairwise interactions have uncovered multistability and chaos Tanaka2011PRL ; Komarov2015PRE ; Bick2016Chaos , few analytical results exist and their collective behavior remains largely unexplored. Here we focus on large systems and obtain an analytical description of the macroscopic dynamics using a partial dimensionality reduction obtained via a variation of the Ott-Antonsen ansatz Ott2008Chaos ; Ott2009Chaos . In particular, the 2-simplex macroscopic dynamics are captured by a combination of two order parameters that capture the degree of synchronization and asymmetry as oscillators organize into two distinct synchronized clusters. Here, clustering refers to multiple phase-locked groups Daido1996PRL ; Skardal2011PRE at different locations on the circle rather than distinct groups of identically synchronized oscillators Pecora2014NatComm . We uncover a novel phenomenon where a continuum of abrupt desynchronization transitions emerge, each at a different critical coupling strength, depending on the asymmetry of the system. Interestingly, no complementary abrupt synchronization transitions occur GomezGardenes2011PRL ; Skardal2014PRE . This continuum stems from an extensive multistability whereby, for sufficiently strong coupling, an infinite number of distinct partially synchronized states are stable in addition to the incoherent state, which is stable for all finite coupling strengths. This multistability indicates the capability of storing a wide array of possible information as different oscillator arrangements. Serving as a minimal model for memory storage, the system captures the critical properties of easily transitioning from an information storage state (i.e., synchronized) to the resting state (i.e., incoherent) Deco2013Trends via abrupt desynchronization. The system then may return to another information storage state with an appropriately chosen perturbation. Moreover, the 1-simplex layer can similarly store information provided that the driving strength from the 2-simplex layer is sufficiently strong. The rich nonlinear dynamics that emerge in this relatively simple extension of pair-wise coupling to three-way coupling highlights the complexity that may arise via simplicial interactions in systems like the human brain and the implications of these behaviors on information storage.

Noting that the dynamics of the 1-simplex layer in Eq. (2) are driven by the dynamics of the 2-simplex layer in Eq. (1), we focus first on the $\theta$-layer our analysis by introducing the generalized order parameters $z_{q}=\frac{1}{N}\sum_{j=1}^{N}e^{qi\theta_{j}}$ for $q=1$ and $2$. Note that $z_{1}$ is the classical Kuramoto order parameter while $z_{2}$ typically measures clustering, which we will see in this system. Using the polar decompositions $z_{q}=r_{q}e^{i\psi_{q}}$ we rewrite Eq. (1) as

[TABLE]

We then consider the continuum limit $N\to\infty$ where the state of the system can be described by a density function $f(\theta,\omega,t)$ which describes the density of oscillator with phase between $\theta$ and $\theta+\delta\theta$ and natural frequency between $\omega$ and $\omega+\delta\omega$ at time $t$. Because the number of oscillators in the system is conserved $f$ must satisfy the continuity equation $0=\partial_{t}f+\partial_{\theta}(f\dot{\theta})$. Moreover, because each oscillator’s natural frequency is fixed and drawn from $g(\omega)$ the density function $f(\theta,\omega,t)$ may be expanded into a Fourier series of the form $f(\theta,\omega,t)=\frac{g(\omega)}{2\pi}\left[1+\sum_{n=1}^{\infty}\hat{f}_{n}(\omega,t)e^{in\theta}+c.c.\right]$, where $\hat{f}_{n}(\omega,t)$ is the $n^{\text{th}}$ Fourier coefficient and $c.c.$ represents the complex conjugate of the previous sum.

We then consider the symmetric and asymmetric parts $f_{s}(\theta,\omega,t)$ and $f_{a}(\theta,\omega,t)$, respectively, of $f(\theta,\omega,t)$ which satisfy $f(\theta,\omega,t)=f_{s}(\theta,\omega,t)+f_{a}(\theta,\omega,t)$ with symmetries $f_{s}(\theta+\pi,\omega,t)=f_{s}(\theta,\omega,t)$ and $f_{a}(\theta+\pi,\omega,t)=-f_{a}(\theta,\omega,t)$. Note that the linearity of the continuity equation implies that if both $f_{s}$ and $f_{a}$ are solutions, then so is $f$. While the asymmetric part $f_{a}$ does not allow for dimensionality reduction, the symmetric part $f_{s}$ does. Noting that the Fourier series of $f_{s}$ is given by the even terms of the Fourier series of $f$, i.e., $f_{s}(\theta,\omega,t)=\frac{g(\omega)}{2\pi}\left[1+\sum_{m=1}^{\infty}\hat{f}_{2m}(\omega,t)e^{2im\theta}+c.c.\right]$, we make the ansatz that each even Fourier coefficient decays geometrically, i.e., $\hat{f}_{2m}(\omega,t)=a^{m}(\omega,t)$ Ott2008Chaos ; Ott2009Chaos . Inserting this and Eq. (2) into the continuity equation, we find that each subspace spanned by even terms $e^{2im\theta}$ collapse onto the same low-dimensional manifold characterized by the condition

[TABLE]

Equation (4) describes the evolution of the complex function $a(\omega,t)$, and thereby $f_{s}$, and depends on the order parameter $z_{1}$. Moreover, $a(\omega,t)$ can be linked to the order parameter $z_{2}$ as follows. First, note that in the limit $N\to\infty$ we have that $z_{2}=\iint f_{s}(\theta,\omega,t)e^{2i\theta}d\theta d\omega$, and after inserting the Fourier series for $f_{s}$ this reduces to $z_{2}=\int g(\omega)a^{*}(\omega,t)d\omega$. To further simplify the relationship we make the assumption that the frequency distribution $g(\omega)$ is Lorentzian with mean $\omega_{0}$ and spread $\Delta$, i.e., $g(\omega)=\Delta/\{\pi\left[(\omega-\omega_{0})^{2}+\Delta^{2}\right]\}$, which has two simple poles in the complex plane at $\omega=\omega_{0}\pm i\Delta$. The integral can then be evaluated using Cauchy’s Residue Theorem Ablowitz2003 by closing the integral contour with a semicircle of infinite radius in the lower-half plane and evaluating at the enclosed pole, yielding $z_{2}=a^{*}(\omega_{0}-i\Delta,t)$. We then evaluate Eq. (4) at $\omega=\omega_{0}-i\Delta$ to obtain

[TABLE]

Next we introduce the rescaled parameters $\tilde{K}=K/\Delta$ and $\tilde{\omega}_{0}=\omega_{0}/\Delta$ with rescaled time $\tilde{t}=\Delta t$, then use polar decompositions, yielding (where we have dropped the $\sim$-notation for convenience)

[TABLE]

Equations (6) and (7) describe the dynamics of the even part $f_{s}$, which falls onto a low dimensional manifold similar to the Ott-Antonsen manifold and describes the evolution of $z_{2}$. However, these equations do not capture the asymmetric part of the dynamics, and moreover they depends on the asymmetric part via $z_{2}$’s dependence on $z_{1}$. As we will see, this reflects the system’s dependence on asymmetry in oscillator arrangements between two clusters.

To close the dynamics we apply a self-consistency analysis to characterize the order parameter $z_{1}$. We first note that, by entering the rotating frame $\theta\mapsto\theta+\omega_{0}t$ we set $\omega_{0}=0$ so that $\dot{\psi}_{1}=\dot{\psi}_{2}=0$. Moreover, by rotating initial conditions $\theta(0)\mapsto\theta(0)+\psi_{1}(0)$ we set $\psi_{1}=\psi_{2}=0$. Equation (2) then implies that oscillators that become phase-locked satisfy $|\omega_{i}|\leq Kr_{1}^{2}$, in which case they relax to one of the two stable fixed points $\theta_{i}=\theta^{*}(\omega_{i})$ or $\theta^{*}(\omega_{i})+\pi$, where $\theta^{*}(\omega)=\arcsin(\omega/Kr_{1}^{2})/2$. These two fixed points correspond to the two clusters that the phase-locked oscillators organize into. Specifically, phase-locked oscillators starting near $\theta=0$ or $\pi$ will end up at the fixed points $\theta^{*}(\omega)$ or $\theta^{*}(\omega)+\pi$, respectively. The phase-locked population is described by the density function

[TABLE]

where the asymmetry parameter $\eta$ describes the fraction of phase-locked oscillators in the $\theta=0$ cluster. On the other hand, oscillators satisfying $|\omega_{i}|>Kr_{1}^{2}$ drift for all time and relax to the stationary distribution

[TABLE]

Next, the order parameter $z_{1}$ is given by the integral $z_{1}=\iint f(\theta,\omega,t)e^{i\theta}d\theta d\omega$, which after inserting the density $f$ as defined by Eqs. (8) and (9) reduces to

[TABLE]

where the contribution from the drifting oscillators vanishes due to the symmetry of $f_{\text{drift}}$. Returning to $r_{2}$, Eq. (6) implies that at steady state we have

[TABLE]

Thus, the macroscopic steady-state is described by Eqs. (10) and (11).

Interpreting these analytical results in the context of numerical simulations allows us to understand novel phenomena that occur in the dynamics of Eq. (1). Beginning with simulations of a system with $N=10^{5}$ oscillators whose natural frequencies are Lorentzian with $\omega_{0}=0$ and $\Delta=1$, we consider initial conditions of varying asymmetry, setting initial phases to $\theta_{i}(0)=0$ and $\pi$ with probabilities $\eta$ and $1-\eta$, respectively. We then begin simulations at $K=16$ and after reaching steady-state slowly decrease $K$ to zero, then restore it slowly to $16$. In Figs. 1(a) and (b) we plot the resulting values for the order parameters $r_{1}$ and $r_{2}$, respectively, for $\eta=1$ (blue, top), $0.95$ (red), $0.9$ (green), $0.85$ (orange), and $0.8$ (purple, bottom). Arrows indicate the direction of $K$ and results corresponding to decreasing and increasing $K$ are plotted as open and filled circles, respectively. As $K$ is initially decreased solutions traverse different partially synchronized states that are determined by $\eta$, until each branch undergoes a discontinuous jump to the incoherent state at different critical coupling strengths in abrupt desynchronization transitions. Importantly, this highlights both a rich extensive multistability (here five different branches are shown, but in the thermodynamic limit an infinite number of such states exist) and a continuum of abrupt desynchronization transitions at different coupling strengths. Partially synchronized branches are characterized by the asymmetry parameter $\eta$, indicating that these complex dynamics arise from different allocations of phase-locked oscillators in the two clusters at $\theta=0$ and $\pi$. Next, as $K$ is restored to its initial value of $16$ no spontaneous transitions back to synchronization occur, indicating no abrupt synchronization to complement the abrupt desynchronization transitions. In Figs. 2(a) and (b) we show that our theory captures these dynamics, plotting in solid curves the theoretical predictions of $r_{1}$ and $r_{2}$, respectively, given by Eqs. (10) and (11). We use the same values of $\eta$ as in Fig. 1 and overlay the results from simulations in open circles, noting excellent agreement.

To better understand the nature of the abrupt desynchronization and multistability phenomena, we investigate the minimal partially synchronized states, i.e., minimal non-zero values $r_{1}^{\text{min}}$ and $r_{2}^{\text{min}}$ allowable for different coupling strengths. We first plot numerical results for $r_{1}^{\text{min}}$ and $r_{2}^{\text{min}}$ in Figs. 2(a) and (b) using black circles. Interestingly, while $r_{1}^{\text{min}}$ appears to decay monotonically with $K$, $r_{2}^{\text{min}}$ remains roughly constant. In fact, the minimal branch $r_{2}^{\text{min}}$ is constant and can be leveraged to analytically describe the branches $r_{1}^{\text{min}}(K)$ and $r_{2}^{\text{min}}(K)$ as well as the corresponding asymmetry value $\eta_{\text{min}}(K)$ and critical coupling strengths $K_{c}(\eta)$ where the abrupt desynchronization occurs. We proceed by inverting Eq. (11), obtaining $Kr_{1}^{2}=2r_{2}/(1-r_{2}^{2})$, which can be inserted into Eq. (10), yielding

[TABLE]

While Eq. (12) appears more complicated than Eq. (10), we note that the coupling strength $K$ has been entirely scaled out of the integral, appearing outside with $(2\eta-1)$. We therefore conclude that if the quantities $\sqrt{K}$ and $2\eta-1$ cancel one another out, i.e., $\sqrt{K}(2\eta-1)$ is constant, it follows that the solution $r_{2}$ in Eq. (12) is independent of $K$. We therefore propose the ansatz $\sqrt{K}(2\eta-1)=\text{const.}$ and use the initial condition $\eta_{\text{min}}(K_{c}(1))=1$, where $K_{c}(1)$ denotes the very first coupling strength where a synchronized state is possible with $\eta=1$, yielding

[TABLE]

Equation. (13) implies that along the minimum branch we have that $\sqrt{K}(2\eta-1)=\sqrt{K_{c}(1)}\approx 2.034$, which can be used in Eq. (12) to compute the minimum branch of $r_{2}$, and in turn $r_{1}$ via Eq. (11), yielding

[TABLE]

In Figs. 2(a) and (b) we plot these minimum branches in solid black curves, which agree with the simulation results. Lastly, by inverting Eq. (13) we find the critical coupling strength $K_{c}$ as a function of $\eta$ where the abrupt desynchronization transition occurs, namely

[TABLE]

In Fig. 3(a) we plot the theoretical prediction of the abrupt desynchronization point $K_{c}(\eta)$ as a solid curve vs observations from direct simulations as black circles, noting excellent agreement.

The dynamics in the 2-simplex layer serves as a minimal model for memory and information storage, capturing a number of critical properties. Each distinct synchronized state corresponds to a specific piece of information, differentiated by the clustering arrangement of the oscillators. Moreover, synchronized state the system can quickly and easily transition to the resting state described by incoherence via the abrupt desynchronized transition. The microscopic properties of this abrupt desynchronization transition are illustrated in Fig. 3(b), where for $\eta=0.9$ the distribution $f(\theta)$ of phases is plotted as $K$ is decreased slowly through the critical value of $K_{c}\approx 6.47$. Before the transition the distribution is asymmetrically clustered about $\theta=0$ and $\pi$ and changes slowly until at $K=K_{c}$ all information is lost as the distribution becomes uniform, representing a resting state.

Lastly, we shift our focus to the $\phi$-layer dynamics, rewriting Eq. (2) as $\dot{\phi}_{i}=\nu_{i}+\kappa\rho_{1}\sin(\varphi_{1}-\theta_{i})+d\sin(\theta_{i}-\phi_{1})$, where $\rho_{1}e^{i\varphi_{1}}=\frac{1}{N}\sum_{j=1}^{N}e^{i\phi_{j}}$. The dynamics of each $\phi_{i}$ can be written in terms of a single sine term, however, due to the fact that all $\theta_{i}$’s are generically different, the target phase and amplitude are all distinct, making self-consistency analyses difficult. However, the microscopic properties may still be understood by comparing $\kappa$ and $d$. Namely, for $\kappa\gg d$ the dynamics are effectively undriven, in which case oscillators that become entrained form a single cluster. On the other hand, for $d\gg\kappa$ the driving term dominates, causing each $\phi_{i}$ to entrain with $\theta_{i}$, thereby mirroring any cluster behavior. This is illustrated in Fig. 4, where we plot $\rho_{1}$ as a function of $d$ given $\kappa=3$ and $K=12$ for asymmetries $\eta=1$, $0.95$, $0.9$, $0.85$, and $0.8$ in blue, red, green, orange, and purple circles, respectively. Dashed horizontal lines denote the $r_{1}$ value for each $\eta$, illustrating that $\rho_{1}\to r_{1}$ in the $d\gg\kappa$ regime and in fact the 1-simplex layer can also store information. However, as $d$ decreases into the $d\ll\kappa$ regime this property is lost as the 1-simplex dynamics no longer reflect the structure in the 2-simplex layer.

The analysis presented here demonstrates how the extension from pair-wise to more general simplicial (specifically, $2$-simplex) interactions in coupled oscillator systems can give rise to a host of complex nonlinear phenomena including information and memory storage. Moreover, these phenomena can be captured and described using analytical methods. In particular, we have characterized a continuum of abrupt desynchronization transitions that occur at different coupling strengths without any abrupt synchronization transitions. This continuum stems from extensive multistability, whereby for sufficiently strong coupling an infinite number of partially synchronized states are stable. These different stable states represent synchronized states organized via different asymmetries into two clusters of entrained oscillators. In addition to highlighting the possible complexity that may arise in coupled oscillator systems with simplicial interactions, we hypothesize that simplicial interactions may give rise to novel nonlinear phenomena in other complex systems.