## 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 (*1*–*3*). 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*, *6*–*8*).

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 oscillators(1)where ϕ_{j} gives the phase of the *j*th 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 (*14*–*16*). 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) distribution(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 by(3)where ϕ_{j} are the phases of the oscillators, *R*_{m} are the phase coherences, and ψ_{m} are the mean phases. Note that *Z*_{0} = 1 follows from the definition of the Daido order parameters. Typically, only the *n* = 1 term is considered and is known as the Kuramoto order parameter. Here, *R*_{1} measures the amplitude of the collective behavior of the oscillator population, with *R*_{1} ≈ 0 indicating desynchrony among the oscillators and *R*_{1} = 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 *R*_{2} > 0 and *R*_{1} ≈ 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 parameters(4A)(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 , 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 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).

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

### Emergence of the scaling

The *m*^{2} 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 such that for *j* = 1, …, *N*. A Taylor series expansion of the Daido order parameters may be written as(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 gives(7A)(7B)which holds whenever the quantity can be considered small and justifies the emergence of the *m*^{2} ansatz we found in both the experimental and simulated data (Fig. 1).

This analysis begs the question of how the COA ansatz and the *m*^{2} 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 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 *m*^{2} 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 *m*^{2} 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.

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 and for the drifting population. Then, the same Taylor series–based argument in Eqs. 6 and 7 considering only the contribution of the locked population gives(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, and 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*)(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 *m*^{2} 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,(10)where η_{i} is a white noise process with 〈η_{i}〉 = 0 and 〈η_{i}(*t*)η_{i}(*t*′)〉2δ(*t* − *t*′)δ_{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 *A*_{ij} = *A*_{ij} = 1(0) if oscillators *i* and *j* are coupled (uncoupled). The degree of the oscillator is then given by . 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 give(11)where *L* is a normalized Laplacian matrix given by *L*_{ij} = δ_{ij} − *A*_{ij}/*d*_{i} and . Our assumptions on the network connectivity dictate that *L* has real eigenvalues that may be ordered λ_{1} = 0 ≤ λ_{2} ≤ … λ_{N} with associated eigenvectors {**v**_{1}, …, **v**_{N}}. 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 *L*^{†}(12)

Allowing for stochastic fluctuations about the deterministic steady state ϕ*, we may compute the expected value as(13)where details of this derivation are given in the Supplementary Materials. If the quantity is small, then our expansion of the Daido order parameters (Eq. 7A) tells us that the *m*^{2} ansatz will provide a good approximation for the phase distribution. Thus, considering Eq. 13, we see that the *m*^{2} ansatz will hold for sufficiently strong coupling strengths for any connected network where 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 *m*^{2} 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.

### 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*)(14A)(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 *m*^{2} ansatz to extract a similar macroscopic model for a large network of noisy, heterogeneous oscillators. In particular, we use the *m*^{2} ansatz as a motivated moment closure to extract a macroscopic model for the order parameter *Z*_{1} 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 (*6*)(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*)(16A)(16B)where *ℑ* denotes the imaginary part of the expression. The Fourier series decomposition of *f* is given by(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 yields(18)where barred quantities are the complex conjugate. In the continuum limit, the Daido order parameters *Z*_{m} are given by(19A)(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 *m*th Daido order parameter, which may be integrated against the frequency distribution *g*(ω) to obtain the composite Daido order parameter *Z*_{m}(*t*) (*22*).

The *m*^{2} 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 *Z*_{m}(*t*) ≈ *Ā*_{m}(*c*, *t*) with *c* = ω_{0} − *i*γ 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 *Z*_{m}(*t*) = *Ā*_{m}(ω_{0} − *i*γ, *t*) allows us to simplify Eq. 18 to an infinite set of coupled equations for the Daido order parameters(20)

Finally, we consider the dynamics of *Z*_{1} and apply our moment closure or , which yields an equation of motion for the Kuramoto order parameter *Z* = *Z*_{1}(21)

Separating the real and imaginary parts gives the macroscopic equations(22A)(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 *m*^{2} 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 *m*^{2} 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 *m*^{2} ansatz provides an accurate description of the macroscopic dynamics (Fig. 4). Moreover, we find that the *m*^{2} ansatz provides a useful upper bound for the collective amplitude *R*_{1}, and the accuracy improves with increased coupling strength. This property may be explained by our result that and that as the entire oscillator population is locked to the mean field.

In the limit of zero noise amplitude (*D*→0), the accuracy of the *m*^{2} 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 *m*^{2} 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 *m*^{2} ansatz should become more accurate. In the next section, we investigate how the *m*^{2} 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 *Z*_{m}(*t*) ≈ *Ā*_{m}(*ω*_{0} − *iγ*, *t*), where ω_{0} − *i*γ 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 *m*^{2} 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 *g*_{c}(ω,γ). Let *h*(ω,γ) = *g*(ω) − *g*_{c}(ω,γ), then the frequency integral (Eq. 19B) becomes(23A)(23B)

The accuracy of the macroscopic model will depend on choosing the dispersion parameter such that the magnitude of the error term is minimized. The *m*^{2} ansatz may then be applied to give the higher-order Daido order parameters with error using the relation .

To compute the error term |*E*_{1}(γ,*t*)|, we recall that *A*_{1}(ω,*t*) may be considered a frequency-dependent version of the Kuramoto order parameter *Z*_{1} (*22*). For oscillators that are entrained to the mean field, we may write(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 *Z*_{1} cancels in the limit of large populations of oscillators (*9*).

For the Kuramoto model, oscillators with |ω| ≤ *KR*_{1} are locked to the mean field with entrainment angle . Therefore, we may rewrite the magnitude of the error integral as(25)

We define as the value of γ such that the error term |*E*_{1}(γ)| is minimized. For *KR*_{1} ≈ 0, we may solve for analytically by expanding *E*_{1}(γ) as a Taylor series in *KR*_{1}. Equating the first-order term to zero gives with *K*_{c} = 2γ_{c}, which matches the results obtained by classical self-consistency arguments (*6*, *9*). Further, numerical solutions of Eq. 25 show that 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, is only weakly dependent on the dynamic variable *R*_{1}(*t*).

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

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

Numerical simulations demonstrate that the macroscopic model (Eq. 26) provides a close approximation to as the coupling strength increases, as shown in Fig. 6 (A and B) for *g*(ω) Gaussian and . More significantly, by plotting the dynamics in the 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 can no longer be considered to be constant in this regime.

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 *m*^{2} ansatz. These two approximations have a certain synergy when applied in tandem, as our method for extracting the dominant frequency mode and the *m*^{2} 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 *KR*_{1} is large, relative to the tails of the frequency distribution, both the *m*^{2} 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 *m*^{2} ansatz, which captures the phase distribution of these systems more accurately.

A simple argument showed the emergence of the *m*^{2} 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 *m*^{2} 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 *m*^{2} ansatz robustly emerged for sufficiently strong coupling strengths. Further, the *m*^{2} 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 *R*^{5} in the collective amplitude equation as compared with the cubic scaling *R*^{3} in the Ott-Antonsen equations (*10*, *17*). We note that a cubic scaling is expected for coupling strengths near the critical coupling strength *K*_{c}, 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 *R*^{5} 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 *m*^{2} 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 *m*^{2} 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 λ = 10^{6} 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 *m*^{2} 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.

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.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).