Research ArticleMATHEMATICS

Macroscopic models for networks of coupled biological oscillators

See allHide authors and affiliations

Science Advances  03 Aug 2018:
Vol. 4, no. 8, e1701047
DOI: 10.1126/sciadv.1701047

Abstract

The study of synchronization of coupled biological oscillators is fundamental to many areas of biology including neuroscience, cardiac dynamics, and circadian rhythms. Mathematical models of these systems may involve hundreds of variables in thousands of individual cells resulting in an extremely high-dimensional description of the system. This often contrasts with the low-dimensional dynamics exhibited on the collective or macroscopic scale for these systems. We introduce a macroscopic reduction for networks of coupled oscillators motivated by an elegant structure we find in experimental measurements of circadian protein expression and several mathematical models for coupled biological oscillators. The observed structure in the collective amplitude of the oscillator population differs from the well-known Ott-Antonsen ansatz, but its emergence can be characterized through a simple argument depending only on general phase-locking behavior in coupled oscillator systems. We further demonstrate its emergence in networks of noisy heterogeneous oscillators with complex network connectivity. Applying this structure, we derive low-dimensional macroscopic models for oscillator population activity. This approach allows for the incorporation of cellular-level experimental data into the macroscopic model whose parameters and variables can then be directly associated with tissue- or organism-level properties, thereby elucidating the core properties driving the collective behavior of the system.

INTRODUCTION

The study of coupled oscillators is important for many biological and physical systems, including neural networks, circadian rhythms, and power grids (13). Mathematical models of these coupled oscillator systems can be extremely high-dimensional, having at least as many degrees of freedom as the number of oscillators as well as additional dimensions for the coupling mechanisms between oscillators. However, the elegant simplicity that emerges at the macroscopic scale in many coupled oscillator populations belies this microscale complexity. Quite generally, these systems demonstrate a phase transition as the coupling between the oscillators is strengthened, leading to the emergence of a self-organized synchronized state (4).

This emergence of a synchronized state from the dynamics of a very high-dimensional dynamical system suggests that a low-dimensional representation of this system should be possible. A major step in this direction was proposed by Winfree in 1967 (5) when he intuitively grasped that for systems of weakly coupled, limit cycle oscillators, the time evolution of each oscillator and the effects of coupling with its neighbors may be described by a single-phase variable. This method, known as phase reduction, reduces the dimension of the coupled system to the number of constituent oscillators and has been used to analyze diverse coupled oscillator systems (1, 68).

In the following years, Kuramoto formalized the mathematical procedure for phase reduction and used it to derive his now famous model for N coupled heterogeneous oscillatorsEmbedded Image(1)where ϕj gives the phase of the jth oscillator, K is the coupling strength, and ωj gives the natural frequency of the oscillator (6). The natural frequencies of the oscillators are typically assumed to be drawn from a distribution g(ω), which reflects the heterogeneity in the oscillator population. The Kuramoto model captures the essential features of many coupled oscillator systems and has been used to study the phase transition to synchrony in detail (9).

However, many biological systems contain thousands of oscillators, making even the phase model a very high-dimensional representation of the dynamical system. A recent breakthrough occurred when Ott and Antonsen (10) discovered an ansatz that can be applied to a family of Kuramoto-like systems to derive a low-dimensional model for the macroscopic behavior of the coupled population. When the ansatz is applied, the long-time behavior of a system of N→∞ heterogeneous oscillators can accurately be described by two differential equations, one for the mean phase of the coupled oscillators and the other for their collective amplitude (11). Despite the hundreds of recent papers that use the Ott-Antonsen dimension reduction procedure, the authors are not aware of any carefully done experiments to test whether this powerful ansatz holds for biological systems.

Here, we test the applicability of the Ott-Antonsen ansatz using a recent experimental data set collected from neurons in the suprachiasmatic nucleus (SCN), the mammalian circadian pacemaker, and through simulations of several models of coupled biological oscillators (12). We find that a core assumption of the Ott-Antonsen ansatz is not valid in our test systems. However, we find that a different, but related, ansatz more accurately describes the data. Using a simple argument, we demonstrate the validity of our ansatz for a wide class of models. We then apply this ansatz to derive a two-dimensional macroscopic model for the population activity of a system of coupled, heterogeneous noisy oscillators. The generality of our procedure should allow for the derivation of low-dimensional macroscopic models of many coupled oscillator systems, allowing for fundamental insights into the core principles driving many biological phenomena.

RESULTS

The development of the Ott-Antonsen ansatz initiated a revolution in the coupled oscillator literature (13). The impact of their ansatz stems from the fact that the macroscopic equations exactly capture all the long-time attractors of the Kuramoto (Eq. 1) and closely related systems (11). The ability to derive strong analytic results has led to its application to a vast array of application areas (1416). Recently, the Ott-Antonsen procedure was applied directly to the study of circadian rhythms for the first time (17).

While an extremely powerful method, the Ott-Antonsen procedure necessarily suffers from several limitations. First, it may only be applied to systems where the interaction between the oscillators is described by a coupling function with a single harmonic (18). Second, the ansatz is not valid for systems whose oscillators evolve with a stochastic component (19). Each of these limitations could severely restrict its applicability to biological systems: Complex coupling forces between biological oscillators often induce higher harmonic components in the model’s coupling function (20, 21), and biological oscillators are invariably subjected to noise (21).

A further limitation of the Ott-Antonsen procedure is one of practicality rather than a formal mathematical restriction. In its most powerful form, the Ott-Antonsen procedure requires the assumption that the distribution of natural frequencies of the oscillators be a rational function g(ω) = a(ω)/b(ω), which is typically taken to be a Cauchy (Lorentzian) distributionEmbedded Image(2)where ω0 is the median frequency and γ controls the range of heterogeneity in the oscillator population. Making this assumption on the frequency distribution is a crucial step in achieving the dimension reduction to the macroscopic model. For more general frequency distributions, the Ott-Antonsen procedure is still mathematically valid, although it produces an infinite set of integro-ordinary differential equations rather than the two-dimensional ordinary differential equation macroscopic model (22). Let us refer to the Ott-Antonsen reduction procedure with the additional assumption of a Cauchy distribution of frequencies as Cauchy Ott-Antonsen (COA).

The ansatz of the COA procedure takes a particularly simple form when written in terms of the Daido order parameters for the distribution of phases of the coupled oscillators. The Daido order parameters (23, 24) are given byEmbedded Image(3)where ϕj are the phases of the oscillators, Rm are the phase coherences, and ψm are the mean phases. Note that Z0 = 1 follows from the definition of the Daido order parameters. Typically, only the n = 1 term is considered Embedded Image and is known as the Kuramoto order parameter. Here, R1 measures the amplitude of the collective behavior of the oscillator population, with R1 ≈ 0 indicating desynchrony among the oscillators and R1 = 1 indicating perfect synchrony.

The higher Daido order parameters add additional details to the description of the phase distribution properties and are often used in the study and detection of cluster synchrony (23, 25). For example, a phase distribution with R2 > 0 and R1 ≈ 0 indicates the presence of two synchronized clusters in the population. Generally, Daido order parameters emerge in the study of phase oscillator systems with higher harmonics in their coupling interactions and phase sensitivity functions (15, 23, 26).

The COA ansatz is a simple geometric relation between the Daido order parametersEmbedded Image(4A)Embedded Image(4B)When the phase distribution of the oscillators is unimodal and symmetric about its mean phase, we expect the mean phase relation ψm = mψ1 to hold generally. Here, we will focus on the case where the phase distribution is approximately unimodal and symmetric. The prediction that Embedded Image, however, is more subtle, and its accuracy has not been evaluated for biological systems.

To test the COA ansatz, we computed the Daido order parameters for a recently published data set measuring the approximate 24-hour oscillations of protein expression in neurons from whole SCN explants (12). Phases were computed from hourly measurements of protein expression in individual neurons over a week-long period as the neurons resynchronized following the application of a desynchronizing perturbation (see Materials and Methods). We examined this data set for evidence of the COA relation Embedded Image at each time point. We found that the phase coherences did not follow this relation (Fig. 1A). In addition, numerical simulations of several different coupled populations of biological oscillator models also reveal that the COA ansatz does not provide a good representation of the equilibrium phase coherences for these systems (Fig. 1, B to D).

Fig. 1 A low-dimensional structure in the phase distribution of coupled oscillator systems.

(A) Experimental SCN neuron data (see Materials and Methods and the Supplementary Materials) (12). (i and ii) Green dots show the phase coherences computed from hourly measurements of cell protein expression in the SCN neurons. The solid black curve shows the relation Embedded Image, and the dashed curve shows the COA ansatz Embedded Image. Inset plots show the circular mean vector of ψmmψ1 across all observations. Bottom row: Histogram (iii) and the first 10 phase coherences (iv) of the phase distribution computed from the data point indicated by the blue star in the top row compared to the phase distribution satisfying the m2 ansatz (black curves). (B to D) Each figure shows a different model simulation: (B) simulation of coupled heterogeneous repressilator oscillators (32, 44), (C) simulation of coupled heterogeneous Morris-Lecar neural oscillators (45, 46), and (D) simulation of coupled noisy modified Goodwin oscillators (47) (see the Supplementary Materials for model details). In each figure (B to D): (i and ii) histogram of the simulated equilibrium phase distribution computed from model simulations for two different coupling strengths [(i), strong coupling; (ii), weak coupling], compared to the m2 ansatz phase distribution. (iii) The first 10 phase coherences for the simulated equilibrium phase distributions for two coupling strengths (green dots, strong coupling; blue squares, weak coupling) compared to the m2 ansatz relation (solid curves).

Instead, in each of these systems, we found that the relationEmbedded Image(5)captures the properties of the phase distribution. We shall refer to this alternate scaling of the Daido order parameters as the m2 ansatz, although we note that this ansatz is equivalent to the Gaussian ansatz used elsewhere (19, 27).

Emergence of the scaling

The m2 ansatz may be derived under more general assumptions than those required by the Ott-Antonsen procedure. Let us consider a population of N coupled oscillators with an equilibrium phase distribution Embedded Image such that Embedded Image for j = 1, …, N. A Taylor series expansion of the Daido order parameters may be written asEmbedded Image(6)

Making use of our assumption that the equilibrium phase distribution is unimodal and symmetric, we have ψm = mψ1, and without loss of generality, we may set ψ1 = 0. Thus, introducing the notation Embedded Image givesEmbedded Image(7A)Embedded Image(7B)which holds whenever the quantity Embedded Image can be considered small and justifies the emergence of the m2 ansatz we found in both the experimental and simulated data (Fig. 1).

This analysis begs the question of how the COA ansatz Embedded Image and the m2 ansatz can both be true. The root of the discrepancy is in the “fat tails” of the Cauchy distribution for the natural frequencies of the oscillator population. The slow decay of the tails of the Cauchy distribution profile results in a significant fraction of oscillators whose phases are not locked to the mean phase but instead drift relative to the population rhythm. This effect keeps the quantity Embedded Image large for any finite coupling strength. However, for natural frequency distributions with exponential tails (for example, Gaussian), the fraction of locked oscillators grows quickly as coupling strength increases, and the m2 ansatz emerges for moderate coupling strengths. Figure 2 (A and B) shows the phase coherences for simulations of the Kuramoto system (Eq. 1) with Gaussian and Cauchy distributions of natural frequencies. Thus, we conclude that the m2 ansatz provides a close approximation for systems with natural frequency distributions with exponential tails, while the COA procedure provides an exact relation for systems with a Cauchy distribution of natural frequencies.

Fig. 2 The emergence of the ansatz depends on the frequency distribution of the oscillators.

The Kuramoto model (Eq. 1) with Gaussian and Cauchy distributions for the natural frequencies of the oscillators, g(ω). (A and B) Relation among the Daido order parameters computed from numerical simulations (circles) and predicted (curves) by (A) the m2 ansatz for Gaussian g(ω) and (B) the COA ansatz for Cauchy g(ω) for increasing coupling strength. Colors indicate different coupling strengths with K that is normalized to the critical coupling strength Kc where partially synchronized solutions emerge (6, 9): K/Kc = 1.1 (red), K/Kc = 1.5 (blue), and K/Kc = 3.0 (green). (C) The fraction of oscillators phase-locked to the mean phase p as a function of normalized coupling strength K/Kc for a Cauchy (dashed green curve) and Gaussian (solid black curve) g(ω).

In fact, we may introduce a correction to our ansatz, which takes into account the presence of phase-locked and phase-drifting oscillators in the population. Let p be the fraction of the population whose phases are locked to the mean phase. Then, the Daido order parameters can be expressed as Embedded Image and Embedded Image for the drifting population. Then, the same Taylor series–based argument in Eqs. 6 and 7 considering only the contribution of the locked population givesEmbedded Image(8)which collapses to Eq. 7B as p→1. In addition, this analysis shows that assuming p = 1 gives a lower bound on the Daido order parameter. In particular, Embedded Image and Embedded Image as p→1. For the Kuramoto model (Eq. 1), we may calculate the fraction of phase-locked oscillators as the coupling strength increases p(K) as (6, 9)Embedded Image(9)

For the Kuramoto model with Gaussian or Cauchy distributions of oscillator natural frequencies g(ω), we may solve for p(K) using Eq. 9. The comparatively slow growth of the fraction of locked oscillators as K increases for the Cauchy distribution relative to a Gaussian distribution of natural frequencies is shown in Fig. 2C.

Complex networks and noise

The simplicity of our derivation makes it clear that the m2 ansatz should hold quite generally. In this section, we characterize its emergence for the case of systems with complex network coupling and intrinsically noisy oscillators. To explore this, we consider a model network of N noisy heterogeneous phase oscillators,Embedded Image(10)where ηi is a white noise process with 〈ηi〉 = 0 and 〈ηi(ti(t′)〉2δ(tt′)δij, where δij is the Kronecker delta. Network connectivity is defined by the adjacency matrix A, and we assume an undirected network such that A is symmetric and Aij = Aij = 1(0) if oscillators i and j are coupled (uncoupled). The degree of the oscillator is then given by Embedded Image. Let the coupling function H be a 2π periodic function and we assume that H′(0) > 0. We note that Eq. 10 is quite general and may be derived in many applications from higher-dimensional, limit cycle oscillator network models under the assumption of weak coupling (28).

We consider the case of strong coupling between the oscillators such that ϕj − ϕi ≈ 0 for all oscillator pairs. In this case, we can linearize about the phase-locked state to giveEmbedded Image(11)where L is a normalized Laplacian matrix given by Lij = δijAij/di and Embedded Image. Our assumptions on the network connectivity dictate that L has real eigenvalues that may be ordered λ1 = 0 ≤ λ2 ≤ … λN with associated eigenvectors {v1, …, vN}. For the linear system (Eq. 11) in the absence of noise (D→0), we may solve for the deterministic steady state ϕ* using the Moore-Penrose pseudoinverse of the normalized Laplacian LEmbedded Image(12)

Allowing for stochastic fluctuations about the deterministic steady state ϕ*, we may compute the expected value Embedded Image asEmbedded Image(13)where details of this derivation are given in the Supplementary Materials. If the quantity Embedded Image is small, then our expansion of the Daido order parameters (Eq. 7A) tells us that the m2 ansatz will provide a good approximation for the phase distribution. Thus, considering Eq. 13, we see that the m2 ansatz will hold for sufficiently strong coupling strengths for any connected network where Embedded Image is finite. In addition, Eq. 13 can be used to study how the emergence of the ansatz depends on the network connectivity, noise strength, and the arrangement of the heterogeneous frequencies in the network (29).

These results are confirmed by numerical simulations of Eq. 10 for the noisy, heterogeneous Kuramoto model [where H(θ) = sin(θ)] with different network connectivity topologies (Fig. 3). In particular, we find that the m2 ansatz provides a quality approximation to the Daido order parameters for both Watts-Strogatz small-world (30) and Barabasi-Albert scale-free (31) network topologies. For each network topology, the accuracy of the approximation increases with the strength of the coupling, as predicted by Eq. 13.

Fig. 3 The equilibrium phase distribution of complex network phase oscillators converges to the m2 ansatz as the coupling strength between the oscillators increases.

(A) Barabasi-Albert scale-free network. (B) Watts-Strogatz small-world network. Circles show the results from simulations of networks of N = 1000 coupled oscillators with noise amplitude D = 1 and oscillator frequencies drawn from a Gaussian distribution with σ = 1. Solid curves show Embedded Image. Colors indicate different coupling strengths as in Fig. 2. Additional details of these simulations are given in the Supplementary Materials.

Macroscopic model

A principal strength of the Ott-Antonsen approach is that the dynamics of the Kuramoto model (Eq. 1) for a large system of coupled oscillators can be reduced to the following two-dimensional macroscopic model (10)Embedded Image(14A)Embedded Image(14B)where ω0 is the median frequency of the oscillators and γ is the dispersion parameter of the Cauchy distribution of natural frequencies (Eq. 2). In this section, we apply the m2 ansatz to extract a similar macroscopic model for a large network of noisy, heterogeneous oscillators. In particular, we use the m2 ansatz as a motivated moment closure to extract a macroscopic model for the order parameter Z1 for the noisy heterogeneous Kuramoto equation (Eq. 10). We consider a fully connected network with coupling function H(θ) = sin(θ). Under these conditions, we may write the system using the Kuramoto order parameter Embedded Image (6)Embedded Image(15)

Following the Ott-Antonsen procedure (10), we consider the continuum limit N→∞ of Eq. 15 and find the continuity equation for the phase density function f(ω,ϕ,t)Embedded Image(16A)Embedded Image(16B)where denotes the imaginary part of the expression. The Fourier series decomposition of f is given byEmbedded Image(17)where c.c. stands for the complex conjugate of the expression and g(ω) is the distribution of natural frequencies of the oscillators. Substitution of the Fourier series for f into the continuity equation yieldsEmbedded Image(18)where barred quantities are the complex conjugate. In the continuum limit, the Daido order parameters Zm are given byEmbedded Image(19A)Embedded Image(19B)using the fact that all oscillating terms in the Fourier series for f integrate to zero except for n = m. Within Eq. 19B, we may regard Ām(ω, t) as a frequency-dependent version of the mth Daido order parameter, which may be integrated against the frequency distribution g(ω) to obtain the composite Daido order parameter Zm(t) (22).

The m2 ansatz applies to the composite, frequency-independent, Daido order parameters. Thus, to apply our moment closure, we need to remove this frequency dependence by obtaining a relation of the form Zm(t) ≈ Ām(c, t) with c = ω0iγ being the “dominant frequency mode.” In the special case where g(ω) follows a Cauchy distribution, the dominant frequency mode approximation is exact with ω0 given by the median frequency and γ given by the dispersion parameter of the frequency distribution (10). For simplicity, we begin by considering g(ω) Cauchy and leave the development of a general approximation scheme for the dominant frequency mode to the following section.

Making the substitution Zm(t) = Ām0iγ, t) allows us to simplify Eq. 18 to an infinite set of coupled equations for the Daido order parametersEmbedded Image(20)

Finally, we consider the dynamics of Z1 and apply our moment closure Embedded Image or Embedded Image, which yields an equation of motion for the Kuramoto order parameter Z = Z1Embedded Image(21)

Separating the real and imaginary parts Embedded Image gives the macroscopic equationsEmbedded Image(22A)Embedded Image(22B)

In previous work, Sonnenschein and Schimansky-Geier (19) derived Eq. 22 for the special case of the noisy Kuramoto model assuming homogeneous oscillator frequencies (γ→0) by using an ad hoc Gaussian moment closure on the phase distribution. The Gaussian moment closure follows the m2 ansatz found here. In agreement with our findings, they found that the macroscopic system (Eq. 22) captured the dynamics of the microscopic noisy homogeneous Kuramoto model accurately, particularly at strong coupling strengths.

Here, we find that the m2 ansatz provides an accurate approximation for the macroscopic dynamics of the noisy heterogeneous Kuramoto model. In Fig. 4, we show the predictions of the macroscopic model (Eq. 22) compared to numerical simulations of the microscopic model in the continuum limit found by using the first 50 moments of Eq. 20 (19). In the case of weak to moderate heterogeneity relative to the noise strength s = γ/D ≤ 1, we find that the m2 ansatz provides an accurate description of the macroscopic dynamics (Fig. 4). Moreover, we find that the m2 ansatz provides a useful upper bound for the collective amplitude R1, and the accuracy improves with increased coupling strength. This property may be explained by our result that Embedded Image and that Embedded Image as the entire oscillator population is locked to the mean field.

Fig. 4 The accuracy of the ansatz for the noisy heterogeneous Kuramoto model.

The equilibrium phase coherence R1 as a function of the coupling strength K for the noisy, heterogeneous Kuramoto model (Eq. 15) for different relative levels of heterogeneity (γ) and noise amplitude (D). (A) s = γ/D = 0.05. (B) s = 0.5. (C) s = 1. (D and E) The transient dynamics of R1 for (D) s = 0.05 and (E) s = 1.0 for different coupling strengths: K = 1.2 (magenta), K = 1.5 (red), and K = 3.0 (blue). In all panels, solid curves show the macroscopic model predictions (Eq. 22), and dashed curves show numerical simulations of the microscopic model in the continuum limit (Eq. 20). Parameters were chosen such that critical coupling strength Kc = 1 for the microscopic model. Insets show curves in the rectangular regions.

In the limit of zero noise amplitude (D→0), the accuracy of the m2 ansatz breaks down under the assumption of a Cauchy distribution of oscillator natural frequencies. This is to be expected, as in the zero noise limit of our system, the COA ansatz has been proven to contain all the long-time attractors of the full system (11). In addition, the Ott-Antonsen approach features an invariance property, which guarantees that states that obey the COA ansatz retain this property for all time (10, 11). Since our ansatz is only approximately obeyed, our method does not share these elegant properties with the Ott-Antonsen approach. However, the use of an approximate ansatz does allow our procedure to be applied to stochastic systems, which are not accurately described by the Ott-Antonsen approach (19).

As previously discussed, the breakdown of the m2 ansatz at weak noise amplitudes (D→0) is related to the fat tails of the Cauchy distribution, which cause the fraction of oscillators locked to the mean field to grow slowly as coupling strength increases. If the natural frequency heterogeneity has less density in the tails of the distribution, our analysis predicts that the m2 ansatz should become more accurate. In the next section, we investigate how the m2 ansatz may be used to derive macroscopic models for systems with strong heterogeneity.

Oscillator heterogeneity

In the derivation of the macroscopic model for the noisy Kuramoto system (Eq. 22), the frequency dependence was removed using an approximation of the form Zm(t) ≈ Ām(ω0, t), where ω0iγ is the dominant frequency mode. As applied in the previous section, this approximation can be made exact when g(ω) follows a Cauchy distribution as the integral may be evaluated as a residue in the lower-half complex plane (10). However, our analysis has shown that the m2 ansatz is best applied to frequency distributions with exponential tails; therefore, we generalize the macroscopic model reduction procedure to allow for more general frequency distributions.

For a general symmetric and unimodal distribution of oscillator frequencies g(ω) with a maximum at ω0, we can think of approximating it with a Cauchy distribution gc(ω,γ). Let h(ω,γ) = g(ω) − gc(ω,γ), then the frequency integral (Eq. 19B) becomesEmbedded Image(23A)Embedded Image(23B)

The accuracy of the macroscopic model will depend on choosing the dispersion parameter Embedded Image such that the magnitude of the error term Embedded Image is minimized. The m2 ansatz may then be applied to give the higher-order Daido order parameters with error Embedded Image using the relation Embedded Image.

To compute the error term |E1(γ,t)|, we recall that A1(ω,t) may be considered a frequency-dependent version of the Kuramoto order parameter Z1 (22). For oscillators that are entrained to the mean field, we may writeEmbedded Image(24)where Ω gives the frequency of the mean field, ρ(ω) describes the collective amplitude, and θ(ω) is the entrainment angle for oscillators with natural frequency ω. When oscillators with frequency ω are locked to the mean field, we have ρ(ω) = 1 (22). The collective contribution of the remaining population of unentrained “drifting” oscillators to the order parameter Z1 cancels in the limit of large populations of oscillators (9).

For the Kuramoto model, oscillators with |ω| ≤ KR1 are locked to the mean field with entrainment angle Embedded Image. Therefore, we may rewrite the magnitude of the error integral asEmbedded Image(25)

We define Embedded Image as the value of γ such that the error term |E1(γ)| is minimized. For KR1 ≈ 0, we may solve for Embedded Image analytically by expanding E1(γ) as a Taylor series in KR1. Equating the first-order term to zero gives Embedded Image with Kc = 2γc, which matches the results obtained by classical self-consistency arguments (6, 9). Further, numerical solutions of Eq. 25 show that Embedded Image decreases quickly as a function of KR and asymptotes to zero for frequency distributions with less density in the tails than the Cauchy distribution (Fig. 5A). Thus, for sufficiently strong coupling strengths, Embedded Image is only weakly dependent on the dynamic variable R1(t).

Fig. 5 Determination of optimal frequency modes.

(A) The optimal frequency mode Embedded Image as a function of KR values, as determined by Eq. 25. (B) The optimal frequency mode Embedded Image when R is given by the long-time asymptotic value Embedded Image. Gaussian g(ω) frequency distribution is shown as a solid green curve and Embedded Image distribution is shown as a dashed black curve. Parameters were chosen such that critical coupling strength Kc = 1 for both distributions.

Therefore, for small perturbations about the synchronized solution in systems with a fixed coupling strength K between the oscillators, we may regard Embedded Image as approximately constant in time. We compute Embedded Image as a function of the coupling strength Embedded Image, where Embedded Image gives the stable phase coherence of the synchronized state. Since the curve Embedded Image can be determined numerically by applying a classical self-consistency approach (6, 9), we can solve for the Embedded Image curve numerically using Eq. 25 (Fig. 5B). We note that each of the distributions we tested gave extremely small error values Embedded Image at the optimal frequency mode.

Following the reduction procedure as given in the last section, with Embedded Image, we find the approximate macroscopic model for the heterogeneous Kuramoto modelEmbedded Image(26A)Embedded Image(26B)

Numerical simulations demonstrate that the macroscopic model (Eq. 26) provides a close approximation to Embedded Image as the coupling strength increases, as shown in Fig. 6 (A and B) for g(ω) Gaussian and Embedded Image. More significantly, by plotting the dynamics in the Embedded Image plane resulting from random perturbations off the synchronized solution, we observe that the macroscopic model provides an accurate approximation of the transient recovery dynamics of the high-dimensional phase model (Fig. 6, C and D). For large deviations from the synchronized state, the approximation breaks down as Embedded Image can no longer be considered to be constant in this regime.

Fig. 6 The macroscopic model for heterogeneous oscillators.

(A and B) The equilibrium phase coherence R1 against the coupling strength K for the Kuramoto model for (A) Gaussian g(ω) and (B) Embedded Image distributions of natural frequencies. Exact solutions obtained from classical self-consistency theory (6, 9) are shown as dashed green curves, and the solution according to the m2 ansatz is shown as solid black curves. Insets show curves in the rectangular regions. (C and D) Plot of the dynamics in the phase plane Embedded Image for perturbations about the synchronized state for K = 3. The dashed black curve shows the predicted dynamics by the macroscopic model, and the solid colored curves show the recovery dynamics for random perturbations off the synchronized state in the high-dimensional phase model. Circles indicate the initial conditions for the transient curves. (C) Gaussian g(ω) and (D) Embedded Image distribution.

Our method for extracting low-dimensional macroscopic models for heterogeneous oscillators can be separated into two key steps: the estimation of a dominant frequency mode and moment closure via the m2 ansatz. These two approximations have a certain synergy when applied in tandem, as our method for extracting the dominant frequency mode and the m2 ansatz are each most accurate for systems where a large fraction of the oscillators are phase-locked to the mean field (p ≈ 1). Specifically, in the regime where KR1 is large, relative to the tails of the frequency distribution, both the m2 ansatz and dominant frequency mode estimation reach their peak accuracy. We find that the resulting macroscopic models provide a close approximation to the dynamics of the full phase model in the vicinity of the equilibrium states. We found similar results when we evaluated the accuracy of the macroscopic model reduction for the full coupled heterogeneous repressilator model in the Supplementary Material (fig. S4) (32). Thus, our method may be applied to derive approximate macroscopic models for a wide class of coupled oscillator systems and study the response to perturbations about the steady states.

DISCUSSION

In the past decade, the powerful ansatz discovered by Ott and Antonsen (10) has been used to resolve many open problems in the coupled oscillator literature and has been applied to an increasing number of application areas (14, 16, 17). Here, we provide the first evaluation of the suitability of the Ott-Antonsen reduction procedure for extracting macroscopic models of real biological networks.

Our examination of a recent experimental data set of circadian oscillator activity (12), as well as simulations of several biological oscillator networks, revealed that these systems did not follow the Ott-Antonsen ansatz. Instead, we identified a new relation, the m2 ansatz, which captures the phase distribution of these systems more accurately.

A simple argument showed the emergence of the m2 ansatz for systems of coupled oscillators, which have a high percentage of the oscillators phase-locked to the mean-field oscillation. We found that the m2 ansatz emerged at moderate coupling strengths for oscillator populations whose frequency heterogeneity has exponential tails. In contrast, the Ott-Antonsen ansatz holds at any coupling strength when the frequency heterogeneity has a Cauchy distribution (polynomial tails). For noisy heterogeneous coupled oscillator systems, the m2 ansatz robustly emerged for sufficiently strong coupling strengths. Further, the m2 ansatz may be used as a moment closure to extract a low-dimensional macroscopic model for noisy heterogeneous oscillator networks.

The low-dimensional system we derive differs slightly from the Ott-Antonsen approach as it produces a term of order R5 in the collective amplitude equation as compared with the cubic scaling R3 in the Ott-Antonsen equations (10, 17). We note that a cubic scaling is expected for coupling strengths near the critical coupling strength Kc, as predicted by the normal form for a Hopf bifurcation (33). Therefore, we expect that our ansatz would overestimate the growth of the phase coherence about the critical coupling strength and may not be an appropriate tool for studying the scaling of the order parameter about the critical coupling. However, as we demonstrated, our approach provides a close approximation to the equilibrium phase coherence as the coupling between oscillators is strengthened.

In the case of human circadian rhythms, several results suggest that models for collective amplitude dynamics should include higher-order terms. For example, higher-order terms in the amplitude growth have previously been required to accurately model the collective amplitude dynamics of the human circadian rhythm in response to a desynchronizing light pulse (34). In addition, the R5 term, which appears in our model, predicts that it should be difficult to increase the amplitude of the circadian rhythm by applying light pulses to an equilibrium circadian amplitude. This is in accordance with experimental results that show that light pulses administered during the day do not significantly affect the circadian amplitude (35). Finally, we note that a previous comparison of two phenomenological van der Pol models for human circadian rhythms showed that the model with higher-order terms better explained human circadian amplitude data (36).

A principal strength of both the Ott-Antonsen procedure and our results is that the parameters and variables of the derived macroscopic models have direct physical interpretations. Therefore, the predictions of the models may be compared with experimental data from the cellular, tissue, and whole-organism levels. For example, Lu et al. (17) made use of the COA ansatz to study jet lag resynchronization asymmetry using readily available data on the mean period of circadian oscillator cells (37, 38). Future work could use this formalism to synthesize cellular-level data on the coupling mechanisms (39), network connections (12), and cellular periods (40) of SCN neurons with behavioral circadian abnormalities observed at the whole-organism level.

To conclude, the m2 ansatz allows derivations of macroscopic models for populations of oscillators with more general frequency distributions and phase-locked behavior than required by the COA ansatz. Our analysis of the phase-locked dynamics of neurons in the mammalian circadian pacemaker suggests that other biological oscillator systems may also be better represented by the m2 ansatz.

MATERIALS AND METHODS

Experimental design

The circadian time series shown in Fig. 1A was collected as described by Abel et al. (12), who generously made their data set publicly available. Briefly, the time series was collected from whole SCN mouse explants cultured for 14 days. The expression of the circadian marker PERIOD2::Luciferase was monitored under a microscope, with bioluminescence measurements collected every hour. On day 6 in culture, tetrodotoxin (TTX) was added to the culture to block neuronal signaling and desynchronize the neurons. The TTX solution was washed away and the culture was allowed to resynchronize. For our purposes, we removed the time points when the TTX solution was added to study the phase distribution of the coupled neurons during resynchronization.

Statistical analysis

The raw bioluminescence data were processed following established methods (40). First, the raw bioluminescence data were de-trended by removing the Hodrick-Prescott baseline trend with a large penalty parameter λ = 106 to minimize loss of the oscillatory signal component. The time-dependent protophase of each oscillator was extracted by dimensional embedding with a 6-hour embedding lag (41). Finally, the time-dependent phase was estimated using the protophase to phase transformation as specified in the Data Analysis with Models of Coupled Oscillators MATLAB toolbox (42, 43).

Details for the mathematical models used in Fig. 1 (B to D) are given in the Supplementary Materials. The estimation of the phase distribution for the in silico data was carried out in the same manner, as described for the experimental data. However, due to the large number of time points available in the simulated data, we used the Hilbert transform to estimate the protophase of the oscillators rather than the dimensional embedding.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/8/e1701047/DC1

Data and mathematical model details.

Emergence of the m2 ansatz for complex heterogeneous noisy networks.

Finding the dominant frequency.

Reduction of limit cycle models to macroscopic models.

Fig. S1. The low-dimensional structure in the phase distribution of coupled oscillator systems from the Abel et al. (12) circadian data set.

Fig. S2. The numerically estimated coupling function Γ(ψ) for the interaction between the repressilator oscillators.

Fig. S3. The variation of the angular frequency induced by the variation of the β parameter in the repressilator model.

Fig. S4. Macroscopic model amplitude recovery dynamics predictions against numerical simulations.

References (4850)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank S. Strogatz, J. Myung, J. Abel, A. Stinchcombe, and A. Cochran for useful discussions. Funding: This work was partially supported by NSF DMS-1412119 (K.M.H. and V.B.), Human Frontiers of Science Program Grant RPG 24/2012 (K.H. and D.B.F.), and NSF DMS-1714094 (D.B.F.). Author contributions: K.M.H. derived the results, wrote the simulations, and analyzed the experimental data. V.B., D.B.F., and K.M.H. conceived the project, interpreted the results and collectively designed the study. V.B., D.B.F., and K.M.H. all contributed to the writing of the manuscript. Competing Interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article