SIAM News Blog

Finding Emergent Low-Dimensional Macroscopic Behavior in Large Systems of Many Coupled Dynamical Units

By Edward Ott

Determining and understanding the macroscopic behavior of a large system consisting of many interacting dynamical parts (“units”) is of interest across a broad range of science. Examples include brain dynamics (where the units are individual neurons), the beating of animal hearts regulated by cardiac pacemaker cells, synchronous signaling in large animal populations (e.g., yeast cells or insects, such as certain species of fireflies), Josephson junction circuits, and chemical oscillators. The macroscopic behavior that typically results in such systems is often described as “emergent” because it cannot be understood purely on the basis of the microscopic dynamics of individual units, but rather “emerges” as a result of their interactions. Perhaps the simplest type of emergent behavior is synchronization in systems of many coupled oscillators (in this case, the oscillators are the units) [13]. For example, cardiac pacemaker cells synchronize throughout a healthy heart, resulting in its rhythmic beating. Seriously malfunctioning hearts, on the other hand, can exhibit more complex emergent behavior in the form of turbulent-like propagating waves of activity (including spiral waves).

While the general subject of emergent behavior is vast, I will focus on a particular (but large) class of models that are often used to address and enhance the understanding of various types of emergent behavior. Somewhat surprisingly, an extremely effective analytical technique uncovers the macroscopic dynamics for the large class of problems to be discussed. I aim to describe this technique here.

To begin, we consider a simple paradigmatic model put forth by Yoshiki Kuramoto in 1975 [2], which provides an excellent illustration of the technique we wish to present. The technique is also capable of treating more complex situations. Kuramoto considers the simplest model of an oscillator, in which the state of oscillator \(i\) is given solely by the phase \(\theta_i\) of its oscillation. The oscillator dynamics are given by \(d \theta_i/dt = \omega_i\), where the constant \(\omega_i\) is the “natural frequency” of the oscillator \(i\). The oscillator population \((i=1,2,\ldots,N)\) is heterogeneous, in that each oscillator has a different natural frequency that Kuramoto assumes to be randomly drawn from some unimodal distribution function \(g(\omega)\). We note that oscillator heterogeneity is natural, e.g., in biological situations where no two units are expected to be exactly the same. Next, Kuramoto pairwise couples his oscillator states to arrive at the model

\[{d \theta_i \over dt} = \omega_i + {k \over N} \sum^N_{j=1} \sin
(\theta_j - \theta_i), \tag1 \]

where the second term on the right-hand side represents the coupling between oscillators with a strength given by the coupling constant \(k\). The main questions addressed by this model are whether or not the oscillators synchronize, and if they do, how strong their synchronization is. Before answering these questions, however, Kuramoto’s complex order parameter can prove useful:

\[ R(t) = {1 \over N} \sum^N_{j=1} \exp[i \theta_j(t)] \equiv r(t)e^{i \psi(t)}, \tag2 \]

Figure 1. Attracting fixed point of the N→∞ macroscopic dynamics of the Kuramoto model, (1), r=|R(+∞)| versus k. kc increases with the spread in the natural frequencies of the oscillators. Image credit: Zhixin Lu.
where \(r(t) = |R(t)|\). For large \(N\) and a uniform random distribution of \(\theta_j\) (i.e., no phase synchronization), \(r \cong 0 \, (r \rightarrow 0\) as \(N \rightarrow \infty\) in this case\()\), while for perfect phase synchrony (i.e., all \(\theta_j\) are equal), we have \(r=1\). Thus, \(r\) quantifies the amount of synchrony for the oscillator population. We can now provide Kuramoto’s main result: for \(N \rightarrow \infty\), the conflict between the desynchronizing tendency (inherent in the spread of natural frequencies \(\omega_i\)) and the synchronizing tendency from the coupling between oscillators is resolved by a so-called “dynamical phase transition” (alternatively, a bifurcation of the macroscopic dynamics), as illustrated in Figure 1. The quantity \(r_\infty = \lim_{t \rightarrow \infty} r(t)\) is the amount of synchronization for the macroscopic attractor of the Kuramoto model, which is \(0\) when \(k < k_c\) and increases continuously as \(k\) increases past \(k_c\), approaching complete synchrony \((r_\infty = 1)\) as \(k \rightarrow \infty\).

The order parameter \((2)\), in addition to characterizing the amount of synchrony, also leads to a useful reformulation of the Kuramoto problem. Specifically, noting that \(\sin \theta = (e^{i \theta} - e^{- i \theta})/(2i)\), \((1)\) can be written as

\[{d\theta_i \over dt} = \omega_i + {i \over 2i} [Re^{-i\theta_i} - R^* e^{i \theta_i}]. \tag3 \]

Since we are interested in large numbers of units \((N \gg  1)\), in order to proceed further it is useful—as Kuramoto does—to make the approximation \(N \rightarrow \infty\), in which case we can characterize the oscillator population at time \(t\) by a distribution function \((\theta,\omega,t)\), such that \(\int^{2 \pi}_0 f(\theta,\omega,t)d\theta = g(\omega)\). And since the number of oscillators is conserved, the continuity equation for the oscillator density \(f\) in \(\theta-\omega\) space becomes

\[{\partial f \over \partial t} + {\partial \over \partial \theta} \left\{ \left[ \omega + {k \over 2i} \left(R e^{-i \theta} - R^* e^{+i \theta} \right) \right] f \right\} = 0, \tag4 \]


\[R = \int \int^{2 \pi}_0 f e^{i \theta} d \theta d \omega. \tag5 \]

Here, \((5)\) comes from noting that, according to \((2)\), \(R\) is the average of \(e^{i \theta}\) over the oscillator population, while \((4)\) results from noting that \(d \omega_i/dt = 0\) (the oscillator natural frequencies do not change) and \(d \theta_i/dt\) is given by \((3)\).

We now consider a certain class of problems—generalizations of the Kuramoto problem—for which the states of coupled dynamical units are given by a single angle variable (i.e., a variable that can be taken modulo \(2\pi\)). I wish to point out the following three main results: 

Figure 2. The invariant manifold M in the space of distributions f(θ,ω,t) and an “orbit” of f, starting at a point (the red dot) on M. Image credit: Zhixin Lu.
Considering the \(N \rightarrow \infty\) limit and evolution of the distribution function \(f\), an invariant manifold \(M\) exists in the space of distribution functions [8]. That is, if \(f(\theta,\omega,0)\) is initialized on \(M\), it continues to evolve on \(M\), as shown schematically in Figure 2.

For appropriate \(g(\omega)\), we can explicitly derive the emergent macroscopic behavior for the system dynamics corresponding to \(f\) on \(M\), often in the form of a low-dimensional set of ordinary differential equations (ODEs) [8]. These two points naturally prompt a key question: Given that initial conditions for \(f(\theta,\omega,0)\) need not be on \(M\), are the two aforementioned points useful? The answer is “yes” because of our third main point [9, 10]: If \(g(\omega)\) is analytic, then all solutions for \(f\) are attracted to \(M\) (i.e., \(M\) is an “inertial manifold”).

That is, using an appropriate metric, the distance from \(f\) to \(M\) approaches \(0\) as \(t \rightarrow + \infty\). This means that we can use the result to analyze all the long-term macroscopic behaviors of these systems, including their attractors and bifurcations.

The specification of the invariant manifold \(M\) comes from the ansatz that the distribution function is of the special form

\[ f(\omega,\theta,t) = {g(\omega) \over 2 \pi} \left\{ 1 + \sum^\infty_{n=1} \left[ \alpha^n(\omega,t)e^{in\theta} + \alpha^{*n} (\omega,t) e^{-in \theta} \right] \right\}, \tag6 \]

where \(\alpha^*\) denotes the complex conjugate of \(\alpha\); i.e., the \(n^{{\rm th}}\) Fourier coefficient in an expansion of the \(\theta\)-dependence of \(f\) is the \(n^{{\rm th}}\) power of \(\alpha(\omega,t)\). Substituting \((6)\) into \((4)\), we find that \((6)\) satisfies \((4)\), provided that

\[\partial \alpha(\omega,t) / \partial t = (k/2) \left[ R(t) \alpha^2(\omega, t) - R^*(t)\right] + \: i \omega \alpha(\omega, t) =
0, \tag7 \]

while \((5)\) yields

\[R^*(t) = \int^{+ \infty}_{- \infty} g (\omega) \alpha(\omega , t) d \omega. \tag8 \]

Determining \(\alpha(\omega,t)\) by solving \((7)\) and \((8)\), since the independent variable \(\theta\) is eliminated, reduces the complexity considerably, as compared to the problem of determining \(f(\theta,\omega,t)\) from \((4)\) and \((5)\). But we can make a further significant complexity reduction by considering particular forms for the natural frequency distribution \(g(\omega)\). The simplest such form is the Lorentzian

\[g(\omega) = {\Delta \over \pi} {1 \over (\omega - \omega_0)^2 + \Delta^2} = {1 \over 2 \pi i} \left\{ {1 \over \omega - \omega_0 - i \Delta} - {1 \over \omega - \omega_0 + i \Delta} \right\}. \tag9 \]

From \((7)\) it follows that analytically continuing \(\omega\) into the complex plane yields the results that \(\alpha(\omega,t) \rightarrow 0\) as \(Im(\omega) \rightarrow - \infty\), and that we can assume \(\alpha(\omega, t)\) to be analytic in the lower half \(\omega\)-plane. Thus, we can evaluate the integral \((8)\) by contour integration to obtain

\[R(t) = \alpha^* (\omega_0 - i \Delta, t). \tag{10} \]

Using \((10)\) and setting \(\omega = \omega_0 - i \Delta\) in \((7)\) yields a single simple ODE for the macroscopic dynamics of the Kuramoto system on \(M\), as characterized by the magnitude of the complex order parameter, \(r(t) = |R(t)|\),

\[dr/dt + (k/2) (r^2-1)r + \Delta r = 0. \tag{11}\]

For \(0 < k < 2 \Delta\), this equation has a single fixed point attractor at \(r=0\), which, as \(k\) increases through the critical point \(k_c=2 \Delta\), bifurcates to an unstable fixed point at \(r=0\) and a fixed point attractor at \(r = \sqrt{1 - (2 \Delta /k)}\). This corresponds to the long well-known result for the Kuramoto model depicted in Figure 1.

While the result \((11)\) is interesting, we have so far not discovered anything new about the Kuramoto model’s behavior. However, ansatz \((6)\) and the technique used to obtain \((11)\) also apply to numerous other problems, yielding many significant new results. The following examples illustrate the great diversity of these results:

  • Modeling walker-induced shaking of the Millennium Bridge [1] (a pedestrian bridge across the Thames River in London), which occurred when the bridge first opened in 2000.
  • Studying why the effect of jet lag can be substantially greater for eastward airplane travel than westward airplane travel by the same number of time zones [6].
  • Modeling brain dynamics of macroscopic synchronized behavior for many interacting, excitable, and firing neurons [7, 11]. Ansatz \((6)\) could be particularly valuable in this area, since the number of neurons in the brain is very large, suggesting that macroscopic reduction might be extremely useful.
  • Studying dynamics that result from different types of coupling between dynamical units, including the effects of time delays along links [4], various network topological properties [12], and consideration of spatiotemporal dynamics in which oscillators are distributed in space and coupled to other oscillators within their local region [3, 5].

In conclusion, the problem of uncovering emergent macroscopic behavior of complex systems, while old, has received much attention in recent years due to its widening relevance to various applications. We hope that the result described in this article will contribute positively to modern developments in these fields.

[1] Abdulrehem, M.M., & Ott, E. (2009). Low Dimensional Description of Pedestrian-Induced Oscillation of the Millenium Bridge. Chaos, 19(1), 013129.
[2] Kuramoto, Y. (1975). Self-Entrainment of Populations of Coupled Nonlinear Oscillators. In H. Araki (Ed.), International Symposium on Mathematical Problems in Theoretical Physics (pp. 420-422). (Vol. 39). Lecture Notes in Physics. Berlin, Germany: Springer.
[3] Laing, Carlo (2011). Fronts and Bumps in Spatially Extended Kuramoto Networks. Physica D, 240(24), 1960-1971.
[4] Lee, W.S., Ott, E., & Antonsen, T.M. (2009). Large Coupled Oscillator Systems with Heterogeneous Interaction Delay. Physical Review Letters, 103(4), 044101.
[5] Lee, W.S., Restrepo, J.G., Ott, E., & Antonsen, T.M. (2011). Dynamics of Pattern Formation in Large Systems of Spatially Coupled Oscillators with Finite Response Times. Chaos, 21(2), 023122.
[6] Lu, Z., Klein-Cardeña, K., Lee, S., Girvan, M., & Ott, E. (2016). Resynchronization of Circadian Oscillators and the East-West Asymmetry of Jet-Lag. Chaos, 26(9), 094811.
[7] Luke, T.B., Barreto, E., & So, P. (2013). Complete Classification of the Macroscopic Behavior of Heterogeneous Networks of Theta Neurons. Neural Computation, 25(12), 3207-3224.
[8] Ott, E., & Antonsen, T.M. (2008). Low Dimensional Behavior of Large Systems of Globally Coupled Oscillators. Chaos, 18(3), 037113.
[9] Ott, E., & Antonsen, T.M. (2009). Long Time Evolution of Phase Oscillator Systems. Chaos, 19(2), 023117.
[10] Ott, E., Hunt, B.R., & Antonsen, T.M. (2011). A Note on the Long Time Evolution of Phase Oscillator systems. Chaos, 21(2), 025112.
[11] Pazo, D., & Montbrio, E. (2014).Low Dimensional Dynamics of Populations of Pulse-Coupled Oscillators. Physical Review X, 4(1), 011009.
[12] Skardahl, P.S., Restrepo, J.G., & Ott, E. (2016). Frequency Assortativity Can Induce Chaos in Oscillator Networks. Physical Review E, 91(6), 060902.
[13] Strogatz, S.H. (2003). Sync: How Order Emerges From Chaos In the Universe, Nature, and Daily Life. New York, NY: Hyperion.

Edward Ott is a Distinguished University Professor in the Department of Physics and the Department of Electrical Engineering at the University of Maryland. He is the recipient of the 2017 Jürgen Moser Lecture. This article is based on the associated lecture delivered in May 2017 at the SIAM Conference on Applications of Dynamical Systems, held in Snowbird, Utah.

blog comments powered by Disqus