^{1}, Thomas M. Antonsen

^{1}and Edward Ott

^{1}

### Abstract

We consider the dynamics of many phase oscillators that interact through a coupling network. For a given network connectivity we further consider an ensemble of such systems where, for each ensemble member, the set of oscillator natural frequencies is independently and randomly chosen according to a given distribution function. We then seek a statistical description of the dynamics of this ensemble. Use of this approach allows us to apply the recently developed ansatz of Ott and Antonsen [Chaos **18**, 037113 (2008)] to the marginal distribution of the ensemble of states at each node. This, in turn, results in a reduced set of ordinary differential equations determining these marginal distribution functions. The new set facilitates the analysis of network dynamics in several ways: (i) the time evolution of the reduced system of ensemble equations is much smoother, and thus numerical solutions can be obtained much faster by use of longer time steps; (ii) the new set of equations can be used as a basis for obtaining analytical results; and (iii) for a certain type of network, a reduction to a low dimensional description of the entire network dynamics is possible. We illustrate our approach with numerical experiments on a network version of the classical Kuramoto problem, first with a unimodal frequency distribution, and then with a bimodal distribution. In the latter case, the network dynamics is characterized by bifurcations and hysteresis involving a variety of steady and periodic attractors.

Emergent behavior in large networks of coupled oscillatory dynamical units is an important issue for many situations arising in diverse scientific and technological settings. These include coupled neurons in the brain, biological oscillators regulating circadian rhythm, Josephson junction circuits, laser arrays, and many others. In this paper, we introduce an “ensemble approach” to the investigation of such systems. Here, by an “ensemble approach,” we mean that, rather than focusing on a particular completely specified system, we instead consider the statistics of a large collection of such systems with similar properties. We show that this approach offers potential numerical, conceptional, and analytical advantages for the study of these systems.

This work was supported by a grant from the U.S. Office of Naval Research (Grant No. N00014-07-1-0734).

I. INTRODUCTION

A. Background

B. Numerical methods

II. UNIMODAL FREQUENCY DISTRIBUTION

A. Formulation

B. Bulk order parameter

C. Steady state

D. Maximum eigenvalue approximation

E. Special case: Uniform in-degree

III. BIMODAL FREQUENCY DISTRIBUTION

A. Formulation

B. Uniform in-degree

C. Dynamics

IV. CONCLUSION

## Figures

In-degree distributions for the Erdös-Renyi and scale-free networks used in this paper. The Erdös-Renyi network’s degree distribution (dot-dashed line) is peaked around 100. Past a minimum degree, the scale-free network takes on a degree distribution (dotted line) of the form , as is more clearly seen in the inset, which is the same plot shown on a log-log scale.

In-degree distributions for the Erdös-Renyi and scale-free networks used in this paper. The Erdös-Renyi network’s degree distribution (dot-dashed line) is peaked around 100. Past a minimum degree, the scale-free network takes on a degree distribution (dotted line) of the form , as is more clearly seen in the inset, which is the same plot shown on a log-log scale.

Eigenspectrum plots for the three networks used in this paper: (a) a directed network with uniform in-degree, (b) an undirected Erdös-Renyi network, and (c) an undirected scale-free network (). In all cases, and . Since the Erdos-Renyi and scale-free graphs are undirected, all eigenvalues in those cases are real.

Eigenspectrum plots for the three networks used in this paper: (a) a directed network with uniform in-degree, (b) an undirected Erdös-Renyi network, and (c) an undirected scale-free network (). In all cases, and . Since the Erdos-Renyi and scale-free graphs are undirected, all eigenvalues in those cases are real.

Bulk order parameter *r* vs. time for systems simulated using the theta formulation (solid line) [Eqs. (1) and (2)], as well as our ensemble formulation (dashed line) [Eqs. (15) and (16)], performed on the networks introduced in Sec. I B: (a) uniform in-degree, (b) Erdös-Renyi, and (c) scale-free. Results were generated numerically using a fourth-order Runge-Kutta integration scheme with fixed time step. Each curve represents a single simulation—no curves are averaged. A time step Δ*t* = 0.1 was used for all theta formulation simulations, save for the scale-free, for which Δ*t* = 0.05 was used, while all ensemble formulation simulations used a time step ten times larger than that was used for the corresponding theta formulation simulations. The width of the frequency distribution was set to Δ = 0.1 and the coupling strength to . Theta formulation curves were shifted horizontally to provide best apparent fit.

Bulk order parameter *r* vs. time for systems simulated using the theta formulation (solid line) [Eqs. (1) and (2)], as well as our ensemble formulation (dashed line) [Eqs. (15) and (16)], performed on the networks introduced in Sec. I B: (a) uniform in-degree, (b) Erdös-Renyi, and (c) scale-free. Results were generated numerically using a fourth-order Runge-Kutta integration scheme with fixed time step. Each curve represents a single simulation—no curves are averaged. A time step Δ*t* = 0.1 was used for all theta formulation simulations, save for the scale-free, for which Δ*t* = 0.05 was used, while all ensemble formulation simulations used a time step ten times larger than that was used for the corresponding theta formulation simulations. The width of the frequency distribution was set to Δ = 0.1 and the coupling strength to . Theta formulation curves were shifted horizontally to provide best apparent fit.

Same as Fig. 3 but with .

Same as Fig. 3 but with .

Same as Figs. 3 and 4 but with .

Same as Figs. 3 and 4 but with .

Long-time-averaged values of *r* vs. *k* for systems simulated using the theta formulation (solid squares) [Eqs. (1) and (2)] and our ensemble formulation (open circles) [Eqs. (15) and (16)] as well as ρ calculated from the transcendental equation (solid line) [Eq. (26)]. The data presented corresponds to the three networks described in Sec. I B. Black points correspond to the uniform in-degree (upper) and scale-free (lower) networks, while the grey points correspond to the Erdös-Renyi network. Also shown as dashed lines are the critical coupling value *k* _{ c }, which is approximately the same for all three networks, and the values of ρ_{max} [Eq. (27)] for the three networks. The same integration scheme was used as for Figs. 3–5. Simulations were generally run for 300 time units for the theta formulation simulations, with averaging done over the last 50 time units, while the ensemble formulation simulations were run until they converged (generally between 200 and 500 time units). Selected points were re-run at smaller time step size and longer simulation runtime to ensure validity.

Long-time-averaged values of *r* vs. *k* for systems simulated using the theta formulation (solid squares) [Eqs. (1) and (2)] and our ensemble formulation (open circles) [Eqs. (15) and (16)] as well as ρ calculated from the transcendental equation (solid line) [Eq. (26)]. The data presented corresponds to the three networks described in Sec. I B. Black points correspond to the uniform in-degree (upper) and scale-free (lower) networks, while the grey points correspond to the Erdös-Renyi network. Also shown as dashed lines are the critical coupling value *k* _{ c }, which is approximately the same for all three networks, and the values of ρ_{max} [Eq. (27)] for the three networks. The same integration scheme was used as for Figs. 3–5. Simulations were generally run for 300 time units for the theta formulation simulations, with averaging done over the last 50 time units, while the ensemble formulation simulations were run until they converged (generally between 200 and 500 time units). Selected points were re-run at smaller time step size and longer simulation runtime to ensure validity.

and vs. for two different values of . When , there is no nonzero intersection of the two curves [thus, no nonzero solution to Eq. (28)].

and vs. for two different values of . When , there is no nonzero intersection of the two curves [thus, no nonzero solution to Eq. (28)].

Bulk order parameter *r* vs. *t* for our uniform in-degree network simulated using our ensemble formulation [Eqs. (15) and (16)] (dashed line) plotted with ρ calculated from Eq. (31) (solid line). The normalized *L* _{2} deviation, *D*, from the manifold given by Eq. (29) is also plotted (dotted line) along with the slope given by Eq. (40) (dash-dotted line).

Bulk order parameter *r* vs. *t* for our uniform in-degree network simulated using our ensemble formulation [Eqs. (15) and (16)] (dashed line) plotted with ρ calculated from Eq. (31) (solid line). The normalized *L* _{2} deviation, *D*, from the manifold given by Eq. (29) is also plotted (dotted line) along with the slope given by Eq. (40) (dash-dotted line).

(a) Bulk order parameter *r* plotted vs. time for a theta formulation simulations [Eq. (1)] on our uniform in-degree network using a bimodal distribution (solid line) and for a simulation using Eqs. (50) and (51) (dashed line). A time step of 0.05 was used for both simulations. The parameters of the simulation were *k* = 40, Δ = 0.5, and ω_{0} = 0.2. (b) A parametric polar plot (*a*, ψ) of the same simulations, starting at incoherent initial conditions (*r* ≪ 1).

(a) Bulk order parameter *r* plotted vs. time for a theta formulation simulations [Eq. (1)] on our uniform in-degree network using a bimodal distribution (solid line) and for a simulation using Eqs. (50) and (51) (dashed line). A time step of 0.05 was used for both simulations. The parameters of the simulation were *k* = 40, Δ = 0.5, and ω_{0} = 0.2. (b) A parametric polar plot (*a*, ψ) of the same simulations, starting at incoherent initial conditions (*r* ≪ 1).

Phase diagram in parameter space showing regions corresponding to different attractor types denoted by I (incoherent steady-state attractor at *r* = 0), SS (steady-state attractor with *r* > 0), and LC (limit cycle attractor corresponding to time periodic variation of *r*). Bifurcations of these attractors occur as the region boundaries are crossed.^{34} The dashed horizontal lines at correspond to the scans of parameter shown in Fig. 11.

Phase diagram in parameter space showing regions corresponding to different attractor types denoted by I (incoherent steady-state attractor at *r* = 0), SS (steady-state attractor with *r* > 0), and LC (limit cycle attractor corresponding to time periodic variation of *r*). Bifurcations of these attractors occur as the region boundaries are crossed.^{34} The dashed horizontal lines at correspond to the scans of parameter shown in Fig. 11.

Long-time behavior of *r* vs. for systems simulated using theta formulation [Eqs. (1) and (2)], plotted in black, and for identical systems simulated using our simplified bimodal ensemble formulation [Eqs. (50) and (51)], plotted in grey, for four different values of : (a) , (b) , (c) , and (d) . When the theta formulation attractors are apparently steady state (SS) or incoherent (I) (in both cases some noise is always present), we discard the first 1000 time units of our simulations (a time step size of 0.05 time units was used), time average of the results over the next 1000 time units,^{35} and plot the averages as solid black squares. When the theta formulation results are apparently limit cycles (noisy), the results are plotted as vertical bars indicating the range of *r* values in the oscillation. Similarly, when the reduced formulation gives oscillation, the vertical range of the grey gives the oscillation range of *r*. Vertical dashed lines represent the region boundaries of Fig. 10. The coupling strength *k* was held fixed at *k* = 40. Simulations were performed on the uniform in-degree network introduced in Sec. I B. Whereas the data in previous figures were obtained by starting the simulations in an incoherent state; simulations in this figure were run twice, once from an incoherent state and again from a coherent initial condition (obtained by pre-running the simulations for large *k*).

Long-time behavior of *r* vs. for systems simulated using theta formulation [Eqs. (1) and (2)], plotted in black, and for identical systems simulated using our simplified bimodal ensemble formulation [Eqs. (50) and (51)], plotted in grey, for four different values of : (a) , (b) , (c) , and (d) . When the theta formulation attractors are apparently steady state (SS) or incoherent (I) (in both cases some noise is always present), we discard the first 1000 time units of our simulations (a time step size of 0.05 time units was used), time average of the results over the next 1000 time units,^{35} and plot the averages as solid black squares. When the theta formulation results are apparently limit cycles (noisy), the results are plotted as vertical bars indicating the range of *r* values in the oscillation. Similarly, when the reduced formulation gives oscillation, the vertical range of the grey gives the oscillation range of *r*. Vertical dashed lines represent the region boundaries of Fig. 10. The coupling strength *k* was held fixed at *k* = 40. Simulations were performed on the uniform in-degree network introduced in Sec. I B. Whereas the data in previous figures were obtained by starting the simulations in an incoherent state; simulations in this figure were run twice, once from an incoherent state and again from a coherent initial condition (obtained by pre-running the simulations for large *k*).

(Color) Polar plot of for a variety of initial conditions for and . The solid lines represent simulations performed using our reduced ensemble equations [Eqs. (50) and (51)] and are color-coded to indicate which attractor each simulation ended on (blue for synchronized steady-state, red for incoherent). The locations of each attractor and of a saddle point are marked by grey dots. The regions surrounding each attractor are blown up in (b) (SS), and (c) (I), with orbits from theta formulation simulations [Eq. (1)] shown in black with transients removed.

(Color) Polar plot of for a variety of initial conditions for and . The solid lines represent simulations performed using our reduced ensemble equations [Eqs. (50) and (51)] and are color-coded to indicate which attractor each simulation ended on (blue for synchronized steady-state, red for incoherent). The locations of each attractor and of a saddle point are marked by grey dots. The regions surrounding each attractor are blown up in (b) (SS), and (c) (I), with orbits from theta formulation simulations [Eq. (1)] shown in black with transients removed.

(Color) Polar plot of for a variety of initial conditions for and . The solid lines represent simulations performed using our reduced ensemble equations [Eqs. (50) and (51)] and are color-coded to indicate which attractor each simulation ended on (blue for synchronized steady-state, red for incoherent). The location of each attractor and of the saddle point is marked by a grey dot. (b) A magnification of the region of interest, with points on the orbit of a theta formulation simulation [Eq. (1)] plotted in black, showing the system starting in the incoherent attractor and escaping to the steady-state attractor. (c) plotted vs. time for the same theta formulation simulation plotted in (b).

(Color) Polar plot of for a variety of initial conditions for and . The solid lines represent simulations performed using our reduced ensemble equations [Eqs. (50) and (51)] and are color-coded to indicate which attractor each simulation ended on (blue for synchronized steady-state, red for incoherent). The location of each attractor and of the saddle point is marked by a grey dot. (b) A magnification of the region of interest, with points on the orbit of a theta formulation simulation [Eq. (1)] plotted in black, showing the system starting in the incoherent attractor and escaping to the steady-state attractor. (c) plotted vs. time for the same theta formulation simulation plotted in (b).

Variance of plotted vs. for fixed values of , and (a) , corresponding to an I attractor and (b) , corresponding to an SS attractor. For each value of , three simulations were run using the same undirected uniform in-degree network with a given value of (generated in the same manner as described in Sec. I B), but using different independent random realizations of the oscillator frequencies. A time step of Δ*t* = 0.05 and a simulation run time of 5000 time units (with averaging done over the last 4000 time units) were used throughout.

Variance of plotted vs. for fixed values of , and (a) , corresponding to an I attractor and (b) , corresponding to an SS attractor. For each value of , three simulations were run using the same undirected uniform in-degree network with a given value of (generated in the same manner as described in Sec. I B), but using different independent random realizations of the oscillator frequencies. A time step of Δ*t* = 0.05 and a simulation run time of 5000 time units (with averaging done over the last 4000 time units) were used throughout.

Two of the graphs from Fig. 11 re-plotted with uniform in-degree data replaced with data corresponding to simulations done on the Erdös-Renyi network were introduced in Sec. I B.

Two of the graphs from Fig. 11 re-plotted with uniform in-degree data replaced with data corresponding to simulations done on the Erdös-Renyi network were introduced in Sec. I B.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content