Model reconstruction from temporal data for coupled oscillator networks

Nature and society brim with systems of coupled oscillators, including pacemaker cells in the heart, insulin-secreting cells in the pancreas, networks of neurons in the brain, cardiorespiratory synchronization, fireflies that synchronize their flashing, chemical reactions, Josephson junctions, power grids, metronomes, and applause in human crowds, to name merely a few.^1–10 The dynamics of coupled oscillators in complex networks have been studied extensively. In particular, networks of interacting oscillators governed by the seminal Kuramoto model^11–13 or its variants are known to exhibit a rich variety of behaviors including spontaneous synchronization, phase transitions, and pattern formation.^14,15 The connections in the network can play a pivotal role in determining these dynamics.

When the governing equations and network topology for a system of coupled oscillators are known, it is fairly straightforward to explore the dynamics of the system using numerical methods and, in special cases, analytical methods. Unfortunately, in many systems, the topology of the interaction network and the intrinsic oscillator properties can be difficult or impossible to observe directly. Imagine, for example, experimental results that track neuronal cell network activity or gene regulatory network activity for a large number of oscillators, giving rise to time series data. One might like to infer models and connectivity information from these time series.

It is crucial to distinguish between functional network connections and structural network connections. Functional connectivity refers to the temporal correlation of the oscillators and is often directly observable. For instance, neurons that fire synchronously are functionally connected. In contrast, structural connectivity, also called network topology, refers to the underlying connections present in the network. For instance, neurons that interact via synapses or neurotransmitters are structurally connected. Unfortunately, these structural connections may not be observable without deconstructing the network.

Network reconstruction is an active area of research, both in theoretical work focused on observations from complex physical systems^16–30 and in applied fields where simultaneously-recorded physiological signals are used to reconstruct physiological networks of interaction.^31–33 Our primary contribution is a new algorithm for solving this inverse problem. We use observed time series data consisting of the phases from a network of coupled oscillators to reconstruct both the network topology, i.e., the structural network connections and the intrinsic oscillator dynamics. For the remainder of this paper, when we refer to inferring or reconstructing the model, we mean both of the aforementioned components: network connections and intrinsic oscillator properties.

In Sec. II B, we review the literature most closely related to our proposed model reconstruction technique. This technique addresses certain challenges of existing methods. One such challenge is that existing methods set up the inverse problem as a linear system that involves many unknown parameters. Consequently, a large amount of time series data is necessary. In contrast, our algorithm results in a system of equations that is smaller and is nonlinear. We solve these using optimization tools designed for training artificial neural networks. As a result, we are able to infer the model with a much smaller amount of data. A second challenge, as mentioned in Refs. 18,22,28, and 29, is that it is difficult to infer the model for networks that are synchronized. For such networks, we explore which perturbations of the synchronized state are needed to enable accurate inference of the model.

The mathematical setting of our study is oscillator networks composed of limit-cycle oscillators. These are also known as phase oscillators, meaning they are characterized by a single phase variable defined on the circle. We explore four different choices of oscillator models, all of which fall into the general framework for uniformly coupled phase oscillators developed in Ref. 34: the classic Kuramoto model,¹¹ a Kuramoto-like model with square-wave coupling function, the Kuramoto-Sakaguchi model,³⁵ and a phase-field reduction of weakly-coupled Hodgkin-Huxley oscillators³⁶ (see  Appendix A). We take all networks to be Erdös-Rényi random graphs.³⁷ We attempt to infer the network topology, the intrinsic frequency of each model oscillator, and the so-called coupling function that specifies the influence of one oscillator on others.

To carry out the model reconstruction, we begin with simulated time series for the oscillator phases and estimated phase velocities. We set up an optimization problem involving the mean squared error for the predicted phase velocities generated by a system of nonlinear differential equations, which can be represented as a convolutional neural network. We then use numerical minimization methods to estimate the optimal parameters, thereby inferring the adjacency matrix for network connectivity, the oscillator coupling strength, the oscillator frequencies, and the Fourier coefficients of the coupling function. This approach is effective for a variety of different networks and model parameters. More specifically, we find the following:

• Accurate inference of the coupling network, frequencies, and coupling functions is possible independent of the model when the system does not synchronize or when sufficient perturbations from synchronization are permitted.

• For small data sets, model reconstruction with our method is more accurate than with alternative methods based on solving linear systems of equations.

• Computational methods designed for training artificial neural networks, such as mini-batch gradient descent implemented in TensorFlow,³⁸ can be used to solve the optimization problem associated with model reconstruction.

• Synchronizing networks can be reconstructed using random phase resets or sufficiently large phase perturbations to a small, randomly chosen set of oscillators.

• Synchronizing networks can also be reconstructed using perturbations to a sufficiently large fixed subset of oscillators.

The rest of this paper is organized as follows. In Sec. II, we introduce the Kuramoto model (a standard model for the dynamics of coupled oscillators), set up the inverse problem for reconstruction, and discuss approaches for solving this inverse problem. In Sec. III, we describe a series of experiments used to test the effectiveness of these reconstruction techniques over a wide range of model parameters. In Sec. IV, we present the results of the aforementioned experiments and discuss perturbation strategies for improving the reconstruction when the system synchronizes. Finally, in Sec. V, we summarize our primary findings and discuss extensions to our methodology.

The Kuramoto model is a standard model for the dynamics of coupled oscillators. We consider $N$ oscillators with phases $θk∈[0,2π)$ for $k=1,2,…,N$⁠, each with an intrinsic frequency $ωk$⁠. These oscillators are coupled through an interaction network with adjacency matrix $A$⁠, with the entry $Akj∈{0,1}$ determining whether oscillators $k$ and $j$ are coupled. The dynamics of $θk$ are governed by the equation

Here, $σkj$ represents the strength of the coupling between oscillators $j$ and $k$⁠, and $Γ(θ)$ represents the coupling function. If the coupling function $Γ(θ)$ is continuous and $2π$-periodic, and $Γ′(θ)$ is piecewise continuous, then it can be represented using a uniformly convergent Fourier series

where the coefficients satisfy $an,bn≤M/n2$ for some $M$⁠. Therefore, these coefficients decay to $0$ for large $n$⁠, and this function can be approximated to arbitrary accuracy using a truncated Fourier series

with sufficiently large $m$⁠. Note that within the model, one can set $a0=0$ without loss of generality, using the change of variables $ωk+∑j=1NσkjAkja0→ωk$⁠.

In Kuramoto’s original formulation,¹¹ the oscillators were globally coupled (⁠$Akj=1$⁠) with coupling strengths $σkj=K/N$⁠, where $K$ scales the global coupling strength. He used the coupling function $Γ(θ)=sin⁡(θ)$ and showed that, for $K>KC$ where $KC$ is a critical value related to the width of the distribution of intrinsic frequencies, the oscillators begin to synchronize, achieving identical phase velocities with phases distributed around the population mean. Analogous results were obtained later for systems with periodic coupling functions possessing only odd harmonics^12,34 and arbitrary complex networks.³⁹ 

We investigate a method for inferring the parameters of the Kuramoto model for phase oscillators on random graphs with uniform coupling strengths $σkj=K/N$⁠. In our experiments, we focus on undirected Erdös-Rényi graphs, where each possible edge in the network is present with a probability of $p$⁠. It is straightforward to extend our approach to arbitrary networks with nonuniform coupling strengths, though this would likely require a larger amount of time series data. The details of our method appear in Secs. II B and II C.

In order to test our method, we generate data according to the following procedure:

1. Generate an Erdös-Rényi network for a fixed connection probability $p$⁠, where the nodes represent individual oscillators with natural frequencies sampled from a normal distribution with mean $μ$ and standard deviation $σ$⁠.

2. Use numerical integration to generate a time series for the evolution of this system for $t∈[0,tmax]$⁠, starting from initial phases drawn from a uniform distribution on $[0,2π]$⁠. For numerical integration, we use an explicit Runge-Kutta method of order 5(4) (Ref. 40) or a stochastic Runge-Kutta method of order 2 (Ref. 41).

3. Compute the phases $θk(tn)$ at times $tn=nΔt$ with time step $Δt$ and for $n=0,1,…,T+1$⁠, where $T=tmax/Δt−1$⁠, to obtain observations that are evenly spaced in time.

4. Estimate the phase velocities $vk(tn)=θ˙k(tn)$ for $n=1,2,…,T$⁠, using central differencing in time with Savitsky-Golay filtering.⁴² We use a window length of $5$ with first degree polynomials.

5. Repeat steps (2–4) $Nres$ times with different uniform random initial conditions to obtain sufficient data during the transient evolution of the system.

See Tables I and II for the network parameters and numerical solution parameters used in this procedure. For these default values, one typically finds that some of the oscillators synchronize with a collective rhythm, while others drift relative to this collective rhythm.

The aforementioned process is intended to produce simulated data mimicking that which experimentalists might collect when observing real world networks. Of course, it may not be possible to control the initial phases in an experiment. We, therefore, evaluate the reconstruction methods for varying $Nres$⁠, as well as for cases in which small perturbations are used instead of different initial conditions (see Sec. IV A).

We now formulate a set of equations whose least squares solution can be used to estimate the natural frequencies $ωk$ of each oscillator, the adjacency matrix $Akj$ of the coupling network, the coupling strength $K$⁠, and the coupling function $Γ(θ)$⁠, from observed phases $θk(tn)$ and phase velocities $vk(tn)$ generated using the method outlined above. This work builds on prior work by Shandilya and Timme,¹⁸ where the adjacency matrix for a network of oscillators is estimated from a time series. We summarize their approach below.

Define vectors $θj=[θ1(tj),θ2(tj),…,θN(tj)]T$ and $vj=[θ˙1(tj),θ˙2(tj),…,θ˙N(tj)]T$ consisting of the oscillator phases and phase velocities. Note that here $T$ corresponds to the matrix transpose. When the coupling strength $K$ and coupling function $Γ$ are known, the phase velocity for oscillator $k$⁠,

is linear in the unknown parameters $ωk$ and $Akj$⁠. Therefore, one can infer both the natural frequencies and adjacency matrix by solving a linear system of $N⋅T$ equations, where $T$ is the number of time steps, and each time step $j$ provides $N$ equations of the form

Here, $L(θj)$ is an $N×(N+N2)$ matrix where each row is determined from Eq. (4) for a particular oscillator, and $x$ is an $(N+N2)×1$ vector with the unknown natural frequencies $ω$ and the entries in the adjacency matrix $A$⁠. The number of unknowns in $A$ can be reduced to $N+N(N−1)/2$ if one assumes $Ajj=0$ and $Akj=Ajk$⁠.

Since the number of equations in the linear system is determined by the number of observed time steps $T$ and the number of oscillators $N$⁠, one can ensure that the system is overdetermined by collecting enough observations so that $T>1+N$⁠. One can then estimate the parameters by minimizing a loss function such as the mean squared error,

where $vj$ denotes the observed value and $v^j$ denotes the predicted value of the velocity.

As long as the system remains far from synchronization, increasing the number of observations $T$ will provide additional linearly independent equations to aid with reconstruction. However, as networks synchronize, i.e., $(vj)k=θ˙k→θ˙$⁠, additional observations become (nearly) linearly dependent. This causes the system of equations to become highly ill-conditioned and makes numerical solutions sensitive to noise and rounding errors. As a result, it is often necessary to perturb the system away from equilibrium to ensure that a sufficient number of linearly independent observations can be collected.

Although the approach of Ref. 18 outlined above is effective, it is limited by the requirement that the coupling function be known a priori. This limitation was addressed in the context of pulse-coupled oscillators in Ref. 28 and phase-coupled oscillators in Ref. 29 by representing the unknown functions as a Fourier series so that it can also be estimated by solving an optimization problem. Unfortunately, one challenge remains. If one represents the coupling function as in Eq. (3), then the system of equations is no longer linear as terms of the form $Akjan$ and $Akjbn$ appear. Pikovsky²⁹ circumvents this issue by defining distinct coupling functions

for each pair of oscillators and by setting $Akj=1$⁠. In this formulation, if oscillators $k$ and $j$ are uncoupled, then $akjn=bkjn=0$ for all $n$⁠. This modification preserves the linearity of the system of equations and allows for the description of a more general class of networks with distinct coupling functions for each pair of oscillators. However, it comes at a cost: the number of unknown parameters increases dramatically from $N+N2$ to $N+2mN2$⁠. Therefore, even longer observation times are required. Additionally, as with the method of Ref. 18, if the system synchronizes, there may be numerical difficulties when inferring parameters. Finally, the large number of free parameters makes this model particularly prone to overfitting. As before, one could assume that connections in the network are symmetric, i.e., that $Γkj(θ)=Γjk(θ)$ and that self edges are not included $Γjj(θ)=0$⁠, allowing for a reduction in the number of parameters to $N+mN(N−1)$⁠.

We propose an alternative approach. We use a single coupling function represented by a Fourier series as described in Eq. (3) for all coupling terms and attempt to infer the $2m$ Fourier coefficients $a=[a1,a2,…,am]T$ and $b=[b1,b2,…,bm]T$ along with the intrinsic frequencies $ω$⁠, adjacency matrix $A$⁠, and global coupling strength $K$⁠. This leads to a total $2m+N+N2+1$ inferred parameters or $2m+N+N(N−1)/2+1$ with symmetry constraints and no self-connections. This is, therefore, a modest increase of $2m+1$ parameters over the case considered by Shandilya and Timme¹⁸ without requiring prior knowledge of the coupling function.

As mentioned previously, this leads to a set of nonlinear equations due to the appearance of terms of the form $KAkjan$ and $KAkjbn$ in Eq. (4). Since $K$ always appears in products involving $an$ and $bn$⁠, these parameters are not structurally identifiable. One could set $K=1$ without loss of generality. Instead, we include $K$ as an inferred parameter and introduce penalty terms to the objective function we seek to minimize,

The inclusion of a penalty function ensures the existence of a nondegenerate local minimum. In this objective function, the first term represents the mean squared error. The remaining terms are penalty terms with weights $λΓ=0.0001$⁠, $λA=10−6$⁠, and $λbd=105$⁠. These hyperparameter values were tuned initially to provide satisfactory reconstruction performance and then kept constant in all subsequent experiments (see also Table III). Additional tuning of these parameters for specific networks could further improve the accuracy of the model reconstruction. The term with $λΓ$ introduces $L2$ regularization, which favors smaller estimates of the parameter values $an$ and $bn$⁠. This is useful for combating overfitting and promoting sparse representations and is analogous to using a Bayesian prior centered at 0 for the Fourier coefficients.⁴³ The term with $λA$ penalizes values in the adjacency matrix that are not close to 0 or 1. The last term with $λbd$ penalizes adjacency values that are negative or greater than 1 to ensure that the estimates fall within the desired range of $[0,1]$⁠. We do not penalize $K$ and instead allow it to be as large as necessary to counteract the parameter shrinkage caused by the penalty terms.

Minimizing Eq. (6) poses a number of challenges. First of all, due to the nonlinearity of Eq. (4), this function may not be convex, and there are no guarantees that local optimization methods will converge to a global minimum. Secondly, when the number of observations $T$ is large, this function may be costly to compute, and minimizing the number of function evaluations is paramount.

Fortunately, the theoretical challenge of nonconvexity does not prevent our method from obtaining consistently reliable model reconstruction, as we show in Sec. IV. The computational difficulties can be addressed by using tools designed for artificial neural networks and implemented in TensorFlow to efficiently compute gradients and perform minimization with a variant of mini-batch gradient descent.⁴⁴ 

To elucidate the connection with artificial neural networks, we comment that Eq. (3) can be viewed as a two-dimensional convolutional neural network with a convolution of size $1×1$ and stride one applied to a tensor of phase differences $Δθi,j,k=θk(ti)−θj(ti)$ with $i=1,2,…T$ and $k,j=1,2,…N$⁠. This network consists of an input layer with dimensions $T×N×N$ (the dimensions of the phase difference tensor), a hidden layer with dimensions $T×N×N×2m$ corresponding to each of the $2m$ harmonics of the form $cos⁡nΔθi,j,k$ and $sin⁡nΔθi,j,k$⁠, followed by a $T×N×N$ output layer representing $Γ(Δθi,j,k)=∑n=1m[ancos⁡n(Δθi,j,k)+bnsin⁡n(Δθi,j,k)]$⁠. A traditional fully connected neural network includes connections between all nodes in consecutive layers, each with a distinct weight. In contrast, a convolutional network uses sparse connections and shared weights that allow it to represent transformations that are equivalent to the mathematical operation known as a convolution. In our network, we only include connections between nodes that correspond to the same phase difference $Δθi,j,k$ and use weights that are independent of the indices $i$⁠, $j$⁠, and $k$⁠. The weights in the first hidden layer are fixed at one, and the biases are fixed at zero. In the output layer, each harmonic is assigned a fixed bias of zero and a variable weight $an$ or $bn$ corresponding to the Fourier coefficients. Figure 1 provides a schematic of this transformation as a neural network. Once the coupling terms have been computed, the rest of Eq. (4) and the resulting loss function [Eq. (6)] are straightforward to compute using vectorized operations.

We initialize the inferred parameters as follows: the adjacency matrix $A$ has initial entries drawn from $N(0.5,1/N)$⁠, the frequencies $ωk$ are initialized from $N(0,1/N)$⁠, the coupling strength $K$ is drawn from $N(1,1/N)$⁠, and the Fourier coefficients $a,b$ of the coupling function are initialized at $0$⁠.

TensorFlow maintains a computational graph for these operations, which allows one to automatically compute gradients. One can, therefore, use gradient descent methods to compute the optimal estimates (⁠$a^n$⁠,$b^n$⁠,$ω^k$⁠,$A^kj$⁠,$K^$⁠) for the parameters. We used a mini-batch gradient descent method with a batch size of $100$⁠; i.e., we randomly assign the data to batches of 100 time steps. For each batch, we calculate the gradient of the loss function and update the estimated parameters by taking a small step in a direction determined from the gradient. Once this has been repeated for all batches (one epoch), we pass through the data again for a total of $200$ or more epochs. We determine the number of epochs experimentally by iterating until the mean squared error no longer decreases. Typically, 300 epochs are sufficient. However, for certain sweeps examining the role of $tmax$⁠, $Nres$⁠, and $σ$⁠, the number of epochs required for convergence is as high as 10 000.

By using small random batches of time steps, the gradient can be estimated quickly. This has the added benefit of introducing stochasticity into gradients, which makes the algorithm less susceptible to getting caught in local minima. We use AdamOptimizer, which is a gradient descent method with momentum and an adaptive learning rate,⁴⁵ with default parameter values for optimization in TensorFlow.

There is minor variability in the success of the algorithm due to the random initialization of the inferred parameters and to the batching of the observed data. In order to ensure an accurate network reconstruction, we can retrain the model with new initial values several times and choose the result that has the smallest mean squared error for the velocity predictions on the validation set, which was not considered during the training process. For any individual model network, this validation error is correlated with the accuracy of the inferred parameters. In all examples below, we attempt to reconstruct each network five times before choosing the best reconstruction.

The method described in Sec. II C produces continuous estimates for the parameters (⁠$a^n,b^n,ω^k,A^kj,K^$⁠), which could be used to predict the dynamics of the system under a variety of initial conditions. However, we are interested in evaluating whether these parameter estimates are an accurate reconstruction of the original model. Before these values can be compared to the model parameters used to generate the data, additional postprocessing is needed. We, therefore, redefine our parameter estimates via the transformations (in order) $K^A^kj→A^kj$⁠, $c0+c1Γ^→Γ^$⁠, where $c0$ and $c1$ are selected to minimize $∫02π|Γ(θ)−Γ^(θ)|dθ$⁠, and $ω^k−c0c1N∑j=1NA^kj→ω^k.$ These transformations are necessary to address the aforementioned identifiability issues with $K$ and to allow $Γ$ to have a nonzero mean.

In the original model, the adjacency matrix entries are either zero or one. However, we treat the entries as continuous variables during optimization and then choose a threshold $ϵ$ so that $Akj<ϵ$ is chosen to be 0 and $Akj≥ϵ$ is chosen to be 1. In practice, one could fix $ϵ$ or select $ϵ$ so that the reconstructed model minimizes the mean squared error for the data set. However, we explore a range of threshold values using receiver operating characteristic (ROC) curves, which measure the quality of a classifier independent of any particular threshold (see  Appendix B). We also report the reconstruction error rates using the value of $ϵ$ that yields the largest $F1$ score for the adjacency matrix reconstruction. The $F1$ score is the harmonic mean of precision and recall and has a maximum value of 1, which corresponds to perfect reconstruction (see  Appendix B). These optimal values are robust, and good performance is typically obtained over a wide range of thresholds (see Sec. IV for details).

In order to validate the robustness of our approach, we test the reconstruction method on data simulated with a variety of networks and parameter values. We explore this high-dimensional parameter space by fixing all of the model parameters except for one and then sweeping the remaining parameter over a wide range of values. See Tables I and II for a list of default values as well as the ranges considered. For each set of parameter values, we compute 30 networks with random initial phases and intrinsic frequencies. We then attempt to reconstruct the underlying model for each network and calculate various performance metrics for each reconstruction.

We evaluate the inferred coupling function $Γ^(θ)$ by comparing it to the true coupling function $Γ(θ)$ via the “normalized difference in area” defined as follows:

This represents the area between the true and estimated coupling function curves, weighted by the area under the curve of the true function. This quantity serves as a measure of the error in the reconstruction. A perfect reconstruction would have a normalized difference in area of zero. The initial estimate $Γ^(θ)=0$ corresponds to a value of 1. Therefore, values significantly lower than 1 indicate progress toward a correct reconstruction.

We compare the inferred intrinsic frequencies $ω^k$ to the true intrinsic frequencies $ωk$ using the “mean absolute deviation” defined as follows:

Values near zero indicate accurate reconstructions. We considered alternative metrics such as relative deviations and the correlation between true and inferred frequencies, but these were less informative because they tend to amplify errors when the intrinsic frequencies are close to zero.

We investigate several methods for evaluating the success of adjacency matrix reconstruction. The accuracy can be inspected visually by plotting the true and reconstructed matrices along with the absolute differences in a grid where values range from $0$ (black) to $1$ (white); see Figs. 2(a)–2(c). Since we are primarily interested in discrete adjacency values, we interpret reconstruction as a classification problem and compute three standard evaluation metrics: the $F1$ score, the classification “error rate” for the connections, and the “area under the ROC curve.” See  Appendix B for precise definitions of these metrics.

Our method can successfully reconstruct models with a variety of coupling functions. Figure 3(a) and Table V show the results for four different coupling functions (see  Appendix A for details), three of which come from popular coupled oscillator models: the Kuramoto model, the Kuramoto-Sakaguchi model, and a phase reduction of the Hodgkin-Huxley model, and a fourth, a square wave, which represents a generic discontinuous function with high frequency Fourier harmonics. For the square-wave coupling function, the accuracy of the reconstruction suffers since the reconstruction includes only five Fourier harmonics, and the true coupling function contains an infinite number of harmonics with amplitudes that decay slowly due to the discontinuities. Despite these limitations, the adjacency matrix was still estimated with a high degree of accuracy.

Counterintuitively, as the number of oscillators increases [Fig. 3(b) and Table VI], we find that reconstruction accuracy of the coupling function improves initially (from $N=5$ to $N=20$⁠). The negative slope in the normalized difference in error is statistically significant with a $p$-value of $0.0001$⁠. This can be explained by the observation that although increasing the number of oscillators increases the number of unknown parameters, it also provides a greater number of pairwise phase differences, which can aid in reconstruction of the coupling function.

In Fig. 4 and Table VII, we consider simulations where the standard deviation $σ$ of the oscillator frequencies ranges from $0.01$ to $1$⁠. The normalized difference in area of the coupling function and the mean absolute deviation of the inferred frequencies both confirm that we can infer the network successfully for $σ≤1$⁠. For larger standard deviations, the quality of the reconstruction suffers slightly, since the large frequencies dominate the coupling terms within the phase velocity relationship.

We also explore the impact of noise on the model reconstruction since both observation noise and dynamic noise would be present in experiments. As in Ref. 18, we introduce dynamic noise using stochastic differential equations with Gaussian noise. We also add Gaussian noise to the observed phases after integrating the governing differential equations. In Tables VIII and IX, we demonstrate that our model reconstruction method is robust to moderate amounts of both types of noise. For both types of noise, we use mean $0$ and standard deviations between $0$ and $1$⁠. Figure 5 shows that dynamic noise with a standard deviation less than $10−3$ helps with the reconstruction of both coupling function and frequencies. This can be explained by the fact that a small amount of noise keeps the system from reaching equilibrium, effectively increasing the duration of the transients that provide useful information about the structure of the network.

We carry out additional parameter sweeps with varying network connectivity $p$ (Table X), maximum simulation time $tmax$ (Table XI), and the number of simulation restarts (Table XII) with all other parameters set using the default values given in Tables I–III. In each case, we are consistently able to reconstruct the coupling function, the underlying network, and the intrinsic frequencies over a wide range of model parameters. Tables V–XVII in  Appendix C illustrate averaged numerical results for 30 trials of each of these parameter sweeps in terms of evaluation metrics such as normalized difference in area, mean absolute deviation, error rate, area under the ROC curve, and a range of thresholds that yield $F1$ scores above $90%$ of the largest value.

As Refs. 18 and 29 point out, when a network remains synchronized, i.e., when $θ˙1=θ˙2=⋯=θ˙N$⁠, one cannot infer the model parameters due to the fact that the observed phases no longer provide linearly independent equations. Furthermore, even with precise knowledge of the adjacency matrix, one could not hope to reconstruct the coupling function $Γ(θk−θj)$ without data over a wide range of phase differences $θk−θj$⁠.

We investigate a method for using small perturbations to introduce brief transients into the dynamics to allow for successful model reconstruction. To test this, we use nearly identical oscillators with frequency standard deviation $σ=0.0001$⁠, which causes the system to quickly converge to a synchronized state for almost all initial conditions. We then initialize $θ(0)=0$ so that the system begins from perfect synchrony. Then, at times $ktmax/Npert$ for $k=1,2,…Npert$⁠, we add a phase perturbation to a subset of the oscillators. Each perturbation causes some of the oscillators to briefly become desynchronized. In this way, the total number of observed phases remains constant, but the fraction of those observations that occur during transient dynamics is proportional to the number of perturbations, $Npert$⁠. The observed phases during these transients provide meaningful data about the structure of the network. Figure 6 demonstrates that, as expected, with too few or too small perturbations, model reconstruction is unsuccessful. However, as the number of perturbations $Npert$ increases, the accuracy of the estimated coupling function, adjacency matrix, and the oscillator frequencies improves.

We explore two methods for selecting which oscillators to perturb: fixed subsets, in which the subset of perturbed oscillators is selected at the beginning of the experiment and these same oscillators are perturbed repeatedly, and random subsets, in which a random subset of the system is selected for each perturbation. Although a perturbation to the phase of a particular oscillator does propagate to the phases of neighboring oscillators through the coupling terms, those perturbations decay quickly and are, therefore, most effective for revealing the local structure of the network. As such, although perturbations to a fixed subset might be more practical in a physical experiment, one must typically perturb a larger fraction of the system to obtain comparable performance to that which is obtained with perturbations to random subsets.

We also consider two types of phase perturbations: phase resets, in which selected oscillators have their phases reset to 0, and phase shifts, in which selected oscillators have their phases modified by adding a random shift $ηpert∼N(0,σpert2)$⁠. Phase resets may be more feasible from an experimental perspective but have the drawback of preserving the mutual synchrony of the subset of oscillators that are perturbed. As such, one will typically need to use random subsets in tandem with phase resets in order to be able to resolve the connections between the perturbed oscillators.

In Figs. 6(a)–6(c), we use phase shifts with $σpert=0.01$ to a fixed subset of size 1, i.e., a single oscillator out of $10$⁠. This gives poor reconstruction regardless of the number of perturbations due to the small size of the perturbation. On the other hand, in Figs. 6(d)–6(f), we perturb a fixed subset (3 out of 10 oscillators) with random phase shifts with $σpert=10$⁠. (These phase shifts are virtually indistinguishable from uniform perturbations $X∼U[−π,π]$⁠.) In this case, the perturbations affect a sufficiently large proportion of the oscillators, and performance begins to improve once there are 5 or more perturbations. The results in Figs. 6(g)–6(i) illustrate that one only needs to perturb 1 out of 10 oscillators to obtain similar performance when the oscillators selected for phase shifts are chosen randomly. In Figs. 6(j)–6(l), we show the results using phase resets to random subsets of 3 oscillators. Again, performance begins to improve dramatically once 5 or more perturbations are used. The tables summarizing the results of these perturbation strategies are provided in  Appendix C.

The amount of transient data necessary for reconstruction is also useful in comparing our method to previous approaches. As discussed in Sec. II B, Pikovsky proposes an alternative method where distinct coupling functions are considered for the interaction of each pair of oscillators.²⁹ This ensures that the system of equations is linear; however, it also increases the number of unknown coefficients significantly. Our approach relies on the assumption that the same coupling function $Γ$ is used for all pairs of oscillators. This is a reasonable approximation for physical systems where the physical mechanisms governing the interactions between oscillators are the same.

Figure 7 shows that when the coupling functions are identical, our method provides more accurate reconstructions of both the network topology and the coupling function with smaller amounts of transient data. As the amount of transient data increases, both methods ultimately achieve perfect network reconstruction, while the approach in Ref. 29 ultimately obtains slightly better estimates for the coupling function. Therefore, our approach would be preferred under circumstances where the amount of data is limited and when the coupling functions are nearly identical.

In this paper, we have designed a method for reconstructing models of coupled oscillator networks including both the network connections and the intrinsic oscillator properties. We hope analysts might investigate (non)convexity of our penalty function and properties of minimizers. These issues aside, after testing our method with many different parameters, we conclude that our algorithm can successfully infer the model.

A common challenge that our method and many related ones encounter is that the procedure fails if the system synchronizes too quickly. One remedy is to have a sufficient number of observations during the desynchronized transients. As we demonstrate, one can ensure that these data are available by repeating experiments with multiple initial conditions or by adjusting the model parameters to inhibit synchronization by instituting large frequency variability (⁠$σ$⁠) or dynamic noise. An alternative remedy is to move the system away from synchrony using perturbations, preferably ones that are physically realizable. The introduction of these perturbations provides useful transient data. By perturbing a sufficiently large subset of oscillators with a large enough change, we are able to infer the model accurately. Our hope is that this method will be adopted by experimentalists and used with experimental data to aid with the construction of interpretable models for the dynamics of networks of coupled oscillators.

Although the numerical experiments outlined here involve Erdös-Rényi networks, the method can be applied more generally to a broad class of networks. Preliminary testing suggests that similar results can be obtained for other topologies such as star, small-world, scale-free networks, and clique networks.

It is also straightforward to extend this approach to other models for coupled phase oscillators such as the Winfree model⁴⁶ or even more general oscillator models such as the Stuart-Landau model⁴⁷ in which both phase and amplitude variations are permitted. Indeed, any system where the unknown functions are periodic can be represented using our technique.

For models containing unknown functions that are not periodic such as the Hodgkin-Huxley model,⁴⁸ a Fourier series representation is not possible. In these cases, one could represent the coupling functions using feed-forward neural networks, which are capable of representing continuous functions to arbitrary accuracy using a finite number of parameters.⁴⁹ Given a sufficiently rich data set, one could still use our approach with backpropagation to learn the structure of the neural network approximation for the coupling function.

This material is based upon work supported by the Mathematics Research Communities of the American Mathematical Society, under the National Science Foundation (NSF) grant (No. DMS-1321794). The project was initiated during the Mathematics Research Community (MRC) on Agent-based Modeling in Biological and Social Systems (2018). A follow-up visit was also supported by a collaboration travel grant from the AMS MRC program. M.J.P. was supported by Hillsdale College. M.-V.C. was supported by The Ohio State University President’s Postdoctoral Scholars Program and by the Mathematical Biosciences Institute at The Ohio State University through NSF DMS-1440386. C.M.T. was supported by the National Science Foundation (Grant No. DMS-1813752) and by Williams College. B.X. was supported by the Robert and Sara Lumpkins Endowment for Postdoctoral Fellows in Applied and Computational Math and Statistics at the University of Notre Dame. We are grateful to Henry Adams, Kelsey Houston-Edwards, and Lori Ziegelmeier for contributions during the formative stages of this work.

The coupling functions investigated are described in Table IV. The first three classes of functions were selected due to their use in popular coupled oscillator models. The last was an example of a generic coupling function with higher order harmonics that was tested to verify that reconstruction of the adjacency matrix is possible even with an imperfect approximation to the coupling function.

The $F1$ score is defined as

Here, $precision$ is the fraction of true positives among all inferred positives, while $recall$ is the fraction of true positives among all positives. The $F1$ score is a value between $0$ and $1$⁠; when reporting the error rate for the reconstructed adjacency matrix, we use the threshold $ϵ=ϵmax$ that corresponds to the largest $F1$ score. We also determine the interval of thresholds that yield $F1$ scores above 90% of this largest value. The width of this interval is reported in Tables V–XVII in  Appendix C as “Interval width $(>90%)$⁠.”

The error rate is defined as the percentage of entries in $A$ that are incorrectly classified. We report the error rate for the optimal threshold $ϵ$ at the end of Sec. II D.

The ROC curve, or receiver operating characteristic, is a parametric curve that uses the classification threshold as a parameter and uses the true positive rate and the false positive rate as the variables.⁵⁰ The area under the ROC curve measures the quality of a classifier independent of any particular threshold. A value of $1/2$ is consistent with blind guesses, and a value of $1$ indicates a perfect classifier since it implies the existence of a threshold for which the rate of true positives is $1$ and the rate of false positives is $0$⁠. Our ROC curves were generated using the scikit-learn package in Python.⁵¹ 

Tables V–XVII report the averaged numerical results of the experimental parameter sweeps described in Sec. III and reported in Sec. IV. The performance metrics used are defined in  Appendix B. Values given in the tables represent the mean plus or minus one standard deviation of the corresponding performance metric over the 30 trials run at that parameter value. We note that some of the distributions, such as the error rate (which is nonnegative by definition), are skewed.