{'fraction_non_alphanumeric': 0.03594667792198657, 'fraction_numerical': 0.010305652280960923, 'mean_word_length': 4.733413994757007, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
KuraNet: Systems of Coupled Oscillators that Learn to Synchronize

Networks of coupled oscillators are some of the most studied objects in the theory of dynamical systems. Two important areas of current interest are the study of synchrony in highly disordered systems and the modeling of systems with adaptive network structures. Here, we present a single approach to both of these problems in the form of"KuraNet", a deep-learning-based system of coupled oscillators that can learn to synchronize across a distribution of disordered network conditions. The key feature of the model is the replacement of the traditionally static couplings with a coupling function which can learn optimal interactions within heterogeneous oscillator populations. We apply our approach to the eponymous Kuramoto model and demonstrate how KuraNet can learn data-dependent coupling structures that promote either global or cluster synchrony. For example, we show how KuraNet can be used to empirically explore the conditions of global synchrony in analytically impenetrable models with disordered natural frequencies, external field strengths, and interaction delays. In a sequence of cluster synchrony experiments, we further show how KuraNet can function as a data classifier by synchronizing into coherent assemblies. In all cases, we show how KuraNet can generalize to both new data and new network scales, making it easy to work with small systems and form hypotheses about the thermodynamic limit. Our proposed learning-based approach is broadly applicable to arbitrary dynamical systems with wide-ranging relevance to modeling in physics and systems biology.

A central theme in the study of complex systems is the emergence of order from disorder, and few formalisms capture this spirit better than networks of coupled oscillators. These systems comprise a set of rotors, each of which is influenced by both its own intrinsic dynamic properties and the extrinsic influence of other rotors in the network. When the intrinsic and extrinsic factors in the network are appropriately balanced, synchrony emerges. That is, oscillators assume a common phase and begin to rotate with a common frequency. The ubiquity of synchronization in nature, from the spiking of neurons [1,2] to the onset of seasonal diseases [3], from the entrainment laser arrays [4] to the release of energy during an earthquake [5], makes oscillatory networks a topic of major practical and theoretical interest in the study of nonlinear and complex systems (for a review of coupled oscillators in complex networks, see [6]).
However, real-world networks are rarely static, and nature abounds with network topologies that evolve towards particular dynamic goals. These adaptive networks add a new wrinkle in the study of oscillatory systems: if the modeler can learn the conditions giving rise to synchrony, can networks be trained to do it on their own? Understanding this more general scenario is crucial to understanding adaptive networks in nature, from modeling changes in brain connectivity [7], to the evolution of flocking behavior [8], and even the emergence of consensus in societies [9]. Moreover, a system which has learned to synchronize provides computational evidence about the conditions which give rise to this collective behavior, lending empirical aid to a theoretically challenging problem. This paper proposes a general framework, KuraNet, for learning in oscillatory systems in the spirit of recent work at the intersection between machine learning and dynamical systems [10][11][12]. KuraNet replaces the standard fixed coupling matrix of the oscillator network with a coupling function whose input data comprises random realizations of intrinsic oscillator features and whose output consists of data-dependent couplings. KuraNet's coupling function takes the form of a deep neural network whose parameters are adjusted by gradient descent. After training, the coupling function stores information about the conditions supporting synchrony in general and not just for a particular sample of disordered features. Unlike previous methods for learning the topology of oscillatory networks [13], our system can generalize to both new oscillator features and network sizes.
We focus on the Kuramoto model [14,15], the archetypical system of coupled oscillators, though our method can be equally applied to other types of oscillatory (and non-oscillatory) dynamics. After a brief description of the approach, we apply it to two general types of collective oscillatory behavior: global synchrony and arXiv:2105.02838v1 [nlin.AO] 6 May 2021 cluster synchrony. In the first case, all phases cohere and oscillators spin at a common frequency. We will show how this technique can easily discover couplings which encourage phase coherence and how these learned functions can shed light on the conditions for synchrony in even highly disordered forms of the Kuramoto model. Cluster synchrony, on the other hand, occurs when oscillators split into separate factions which are internally coherent but mutually desynchronized. In contrast, earlier work [16][17][18] on cluster synchrony has focused on the related but distinct case in which clusters behave independently of each other instead of maintaining strict inter-group separation. We show how the approach can be used to split populations into separate clusters based on class labels and how such a method can be used for the segmentation of high-resolution images. We will conclude with a summary of key theoretical and practical problems to which the framework might be applied in the future.

I. BACKGROUND
An oscillator is a dynamical system with a stable, periodic limit cycle. Malkin [19], Winfree [20], and others showed that, under weak coupling conditions, n oscillators, each associated to a node on the graph G = (N, E), can be reduced to a system of phases, θ(t) ∈ T n , on the n-torus. The most famous such "phase-reduced" system of oscillators is the Kuramoto model, which we consider in the following form: The dynamics of an oscillator, i, are determined by both intrinsic and extrinsic factors. Here, the intrinsic factors comprise the node features: natural frequencies, ω i , external field strengths, b i , and transmission delays, τ i . The external influence on an oscillator is determined by interactions with other phases in the rest of the network, θ j for j = i, via the matrix of weights or couplings, K ∈ R n×n associated to G.
Disorder is introduced by placing probability distributions, f, g and h, on ω i , τ i and b i , respectively. We write the joint distribution as D and denote its sample space, the node feature space, by X . Global synchrony in the system is typically measured by the so-called global order parameter, so that R ≈ 0 when the oscillators are out of phase and R = 1 when oscillators are in phase (i here denotes the imaginary unit, not an index). Below we will instead work with the circular variance, V = 1 − R.
The case of g = h = δ(0) = f with f symmetric and unimodal has been extensively studied, especially for mean field couplings, K ij = /n, for > 0 (see [21,22] for excellent reviews). In this paradigmatic case, there exists a critical coupling, c , such that V → 0 only when > c . The critical coupling is inversely proportional to the variance of f so that stronger coupling is required for more heterogeneous intrinsic frequencies. The case of more general disorder, the focus of the current study, has received less attention and disordered models are known to produce esoteric and poorly understood behaviors, like quasi-glassiness [23,24], hysteresis [25] and multistability [26].

II. KURANET
Various approaches to optimizing Kuramoto network topology have been proposed [13,27,28], though these techniques generally involve the optimization of a single coupling matrix for a given realization of random node features. Consequently, anything learned for this single instance does not necessarily generalize to networks with new disordered features or to networks of a different size. Such techniques are also limited in their ability to model adaptive oscillatory systems in the real world, which must presumably learn to synchronize in the presence of variable forms of disorder and not just a single example. Moreover, the pioneering approaches of [13] and others typically rely on greedy, stochastic methods, whereby single network edges are turned on or off after each of many repetitions of the dynamics, making optimization prohibitively slow for large systems.
Inspired by the growing body of research at the intersection of dynamical systems and machine learning [11,12,29], we propose a new technique called Ku-raNet for building oscillatory systems which can learn to synchronize ( Figure 1). The key difference between our approach and earlier techniques lies in our replacement of the traditional coupling matrix, K, with a differentiable coupling function, K µ : X 2 → R. This function maps from pairs of intrinsic node features to the coupling weight between those nodes. The coupling function takes the form of a deep neural network whose learnable parameters, µ, are optimized to store knowledge about which couplings are appropriate for any realization of random node features. As a result, the coupling function can easily be transferred to new data. Moreover, since K µ is applied to all pairs of nodes in parallel, the system can be automatically scaled up to much larger networks than those observed during training.
KuraNet is trained by stochastic gradient descent on random batches of node features, X ∈ X n (Figure 1A). All n(n − 1) distinct pairs, (x (i) , x (j) ), from a batch are passed in parallel through the coupling function, K µ : X 2 → R ( Figure 1B). The result is a datadependent coupling matrix with elements K µ (x (i) , x (j) ) which then participates in a Kuramoto dynamics ( . First, a batch, X ∈ R n×m , of random intrinsic oscillator features is sampled from the node feature space, X . Here, we show a batch of n = 10 oscillators each with feature, x (i) ∈ R 2 , distributed according to a two-component gaussian mixture. Two oscillator features, x (1) and x (2) , are shaded in pink and followed throughout KuraNet processing. (B). Then, all pairs (x (i) , x (j) ) are passed in parallel (only one pair shown) through the coupling function, Kµ, which takes the form of a deep neural network whose parameters, µ, comprise synaptic weights and biases. (C). The output of the function on the pair of samples (x (i) , x (j) ) becomes the coupling strength Kµ(x (i) , x (j) ) in a Kuramoto model, each of whose i = 1, . . . , n oscillators inherits the intrinsic feature x (i) . A dynamics is simulated for the resulting Kuramoto model and evaluated according to a loss function whose output is used to adjust the parameters, µ, of the coupling function.
. Given an initial condition, θ 0 , Eq. 1 is then solved with any off-the-shelf numerical method, and the solution θ(t; θ 0 , K µ ) is evaluated by where : T n → R is a differentiable instantaneous loss function, t 0 is a burn-in time, and K µ (X) is shorthand for the data-dependent matrix of couplings. It will sometimes be desirable to constrain this matrix to lie in some feasible set by restricting µ ∈ M. The coupling function, the equation of motion for θ and the instantaneous loss function, , are all differentiable in µ so we may approximate by gradient descent, where D n denotes the joint distribution over the batch of node features.
Depending on which type of dynamics we wish to encourage, we will employ two different loss functions. To encourage global synchrony (see Section Results III A), we can simply minimize the circular variance, (θ) = V (θ). To encourage cluster synchrony, we must first assume there exists a function Y : X → {1, . . . , k}, mapping from the node feature space to a set of discrete labels. This labeling associates to the i th oscillator a target cluster via the label of its feature, l = Y (x (i) ). Our goal is then to ensure the phases of all oscillators within a cluster are synchronized and all phases between clusters are maximally desynchronized. We therefore use an instantaneous loss function of the form where V l is the circular variance of the l th group, and S is another loss designed to encourage splayness [4] among target groups in the sense that the mean phases within the groups equidistribute the unit circle. Let θ l be the average phase of oscillators having features x j such that Y (x j ) = l . Then, we set which is minimized exactly when { θ l } k l=1 are splay [30]. The sum in Eq. 6 is a generalization of the global order parameter, R, to higher circular moments. These moments pluralize the notion of synchrony so that phases are considered synchronous and maximize the g th moment as long as they all point in the same direction modulo 2π/g. Conversely, the sum is minimized when the phases are out of phase modulo 2π/g; i.e., they are splay. More details can be found in [30,31]. Note that all losses reported below are all normalized to be between 0 and 1. For the cluster synchrony experiments, the node feature space, X , only controls the labels of oscillators and we assume that f = g = h = δ(0), so that the oscillators are equivalent to XY spins [21, Section IV.B.1].

A. Global synchrony
Here, we learn coupling functions, K µ , which minimize the circular variance, V , of Kuramoto phases for three standard forms of Eq. 1 depending on which parameters are disordered: 5,15), where U and N B are the uniform and negative binomial distributions, respectively. The discrete samples from the negative binomial distribution were used to index the timesteps of an Euler update used to solve Eq. 1. Note that Case III uses a decentered natural frequency distribution since the presence of time delays introduces a nontrivial dependence between system behavior and the natural frequency mean [32]. In all cases, we forced K µ (X) to be symmetric and positive.
To make the optimization non-trivial, we follow earlier work by [13] by restricting the underlying weighted graph associated with K µ (X) to have a fixed, small average degree, in this case d = 1 n ij K µ (x (i) , x (j) ) = 1. This forces KuraNet to be judicious about its limited coupling resources: only those oscillators whose intrinsic parameters are crucial for synchrony deserve to be strongly coupled.
To emphasize KuraNet's ability to generalize across both data and scale, we trained the system on relatively small networks (n = 100) and then generalized to both new data and larger network sizes (n = 10, 000). This gives rise to three testing regimes: Data Testing, where the model is exposed to new node features from the testing set for networks of size n = 100; Size Testing, where the system is exposed to a single network of size n = 10, 000 comprising the whole training set; and Data+Size Testing, where the model is exposed to a single network of size n = 10, 000 comprising the whole testing set. These testing regimes measure the degree to which KuraNet can abstract beyond the observed training data and beyond finite-size effects. Each case was optimized for 10 random seeds and the best values are shown here. Further training details are found in the SI Methods section.
The circular variance, V , achieved in various testing regimes are collected in Table I. Parenthetical values show the corresponding results for an untrained control system, which nevertheless has the same average degree d as the trained one. For all three cases and in all three testing regimes, the optimized model performs substantially better than the control model. The most striking results are achieved in Case I, where the optimized system achieves nearly perfect synchrony in all three testing regimes (Figure 2 A).
To help explain the results of our optimization, we introduce a generalized version of the graph-induced autocorrelation function proposed by [13], where x = 1 n n i=1 x (i) ∈ X . Eq. 7 measures the correlation between random node features with respect to the couplings, K µ (x (i) , x (j) ). Plotting V versus C x for each of the three optimizations reveals that graphinduced anti-correlation between nodes is strongly associated with synchrony (linear at r 2 = 0.81, r 2 = 0.64) in Case I and Case II, but not Case III (r 2 = 0.14). Hence, synchrony in the former cases is promoted when the limited coupling resource is spent between nodes with disparate features, effectively reducing the heterogeneity of the system. A similar result was found by [13], but only for the equivalent of our Case I. Case III is evidently rather different, but we can gain some insight by plotting the learned coupling function. Discretizing the sample space of ω and τ into 30 bins and plotting the average fan-in weight k j = 1 n i K µ (x (i) , x (j) ) associated to each bin over 100 batches shows a marked peak at (ω * , τ * ) = (2.00, 21.00). Notably, the peak occurs at the lowest intrinsic frequency and at a time delay slightly to the right of the mode of h (τ = 15). This was consistent for all of the best random seeds for this case. The mean of the intrinsic frequency distribution sets an effective time-scale for the network and it is known that delays which are long compared to this time-scale tend to impede synchrony [32]. Our result suggests that this impediment can be overcome when weights simultaneously emphasize the longest time-scale (i.e. lowest intrinsic frequencies) and those delays which are as long as possible without being too far from the mode of h. We simply take the results from Case III as a computational hypothesis about disordered time delays in the Kuramoto model, which we leave for future investigation.

B. Cluster synchrony
Next, we turn to the case of cluster synchrony, in which oscillators split into internally coherent factions. Cluster synchrony has attracted a great deal of recent interest and has been observed, for example, in flocking animals [8], models of visual processing [33], and power grids [34]. Here, we are interested in the stronger condition in which oscillators not only synchronize into clusters but also where these clusters are maximally separated in the sense that the cluster phase averages, θ l , are equally spread out across the unit circle. We investigated optimized cluster synchrony for two types of data: simple mixture distributions on the plane, and images taken from the Berkeley Segmentation Data Set (BSDS) [35].
For simple mixture distributions, we used several standard distributions from Python's sklearn software: two interlinked crescent moons ("Moons"), two concentric circles ("Circles"), two interlinked spirals ("Spirals"), as well as nine Gaussian mixture models (GMMk) with between k = 2, . . . , 10 means spaced evenly around a ring of radius 10 (results here only show k = 5, but more are available in the SI Detailed Clustering Results). Node features were the coordinates (p, q) of the corresponding datum and labels came from the labels of mixture components.
As before, we are interested in both generalization to new data and to larger networks, resulting in the three testing regimes described above. Mixture distributions used training and testing network sizes of n = 100 and n = 10, 000 as before. The coupling function architecture was identical to that of Section III B, except that we now permit negative couplings. The coupling matrix was still normalized so that each unit had an (absolute) in-degree of 1.0 before being symmetrized by K ← .5(K + K T ). For details, see SI Methods.
After training, we found that the model was able to generalize to correct splay dynamics across all testing regimes far better than an untrained model with random parameters but with the same degree normalization (Table II). A fuller picture of the trained system is provided by plotting the limiting phase of a single oscillator (five-class Gaussian mixture shown in Figure 3A) as its feature vector is varied across a grid of values. For two cluster data (Moons, Circles, Spirals), the randomized model tended to synchronize globally for a loss of .5, while the optimized model achieves a loss of nearly 0. The optimized system also achieved nearly perfect cluster synchrony in the multi-class cases (e.g. GMM5), where the random model's partial synchrony only produced a loss of about an order of magnitude worse. Videos of the dynamics can be found in the SI.
Plotting couplings sorted by cluster index for a sample of 1,000 training oscillators revealed that the system had learned in all cases a distinctive mixture of block and diffusive connectivity ( Figure 3B). Couplings between oscillators in different groups were uniform as a function of the group indices, and these couplings grew weaker and more negative as the group features grew more distant The phase landscape in X = R 2 is plotted by recording the asymptotic phase of a single oscillator, i, in a fixed batch (small circles) as its 2-D feature was varied over a grid of 100 x 100 values. (B) The trained system learns to produce coupling matrices with a clear block-diffusive structure whose cells correspond to the true labeling function on X . (C) However, not just any block coupling structure can produce this effect, since the optimized model (right bars) achieves a lower total loss than both a randomized model (left bars) and a block control model (middle bars) with positive intra-cluster couplings and negative inter-cluster couplings. (D,E,F) The differences among these models is all the clearer when plotting terminal phases (small dots) and cluster mean fields (colored arrows). Here, phases and mean fields are colored according to ground truth clusters, not angle. While the control model can easily exhibit strong intra-group synchrony (E, only the balance of positive and negative weights in the learned diffusive couplings of the optimized model (F) can also produce reliable splay configurations.
in the plane. Earlier work [16][17][18] has implicated one of these properties, block coupling, in the emergence of cluster synchrony, although these studies did not include our extra condition of mutual desynchrony among clusters. To compare our optimized couplings to coupling strategies in this earlier work, we examined an additional control model with block couplings normalized exactly as in the optimized case but with a fixed intra-group (inter-group) coupling of b > 0 (b < 0). The value of b was chosen to such that the optimized and control settings would have the same maximum absolute coupling strength, max i,j |K ij |, so that attraction and repulsion between groups would be on roughly the same scale in both cases. We found the block control could successfully synchronize each cluster far better than the random setting ( Figure 3C), although the independence of clusters prevented any rigid splay behavior, particularly for k > 2. This is especially evident when the mean fields ( Figure 3 D-F) are plotted for the true clusters in each condition. This result suggests that the learned diffusivity in the optimized system enforces cluster separation by tuning the balance between the relatively few intra-group interactions and relatively many inter-group interactions which prevail in the case of k > 2.
For the image experiments, we chose eight images from BSDS. Here, KuraNet functions as a form of image segmentation (Fig. 4). Various other oscillator-based models of image segmentation have been proposed [36,37], though, to our knowledge, the current work is the first end-to-end differentiable model trainable by gradient descent. Similar to the mixture distribution case, node features include the pixel location location (p, q), although we now also include the the RGB channel intensities I r , I g , I b ∈ [0, 255] at the location (p, q). Labels came from the ground truth segments of the BSDS data set. Training was identical to the mixture distribution case except that training and test sets were considerably larger owing to the image sizes (∼ 150, 000 pixels). Note that KuraNet learns to segment one image at a time (as opposed to multiple images). We speculate on ways to achieve this larger goal in Section IV.
Performance on images was favorable (Table II, last  row; full results, SI Detailed Clustering Results), and Fig.  4 shows the results on the full images (training plus test sets). Note that the images in Figure 4 represent the dynamics on very large networks (again, ∼ 150, 000 nodes) even though KuraNet was only ever exposed to dynamics on networks three orders of magnitude smaller. Discretizing the limiting phase configuration (Fig. 4, second column) using k-means clustering (Fig. 4, third column) and measuring performance using the standard symmetric best dice score 0 ≤ SBD ≤ 1 [38] shows that KuraNet achieves an accuracy of .7001 vs a value of .3719 (slightly above chance level) attained by k-means on raw pixel in-  Table II). (First column) Original image. (Second column) Limiting phase of KuraNet on these images after training (hue represents phase). (Third column) k -means clustering is used to discretize the phase into segments. (Fourth column) Corresponding k-means segmentations using raw pixels as data instead of KuraNetproduced phases. (Fifth column) Ground truth segmentation. Note that the correspondence between predicted segments and ground truth segments does not rely on the absolute labeling of either segmentation (here, the color), but rather the existence of a bijection between the two sets of labels which preserves the relationship between any two pixels' labels (i.e., same vs. different) (see SI Methods.) tensities (fourth column). These results cast KuraNet as an interesting new direction for image processing using complex systems.

IV. DISCUSSION
Recent work at the intersection of machine learning and dynamical systems has typically taken a somewhat model-free approach, estimating couplings [29] or learn-ing dynamics [10] while assuming little about the underlying equation of motion. Here, we demonstrate the value of the model-based setting, whereby the parameter space of a given dynamics of interest, Kuramoto flows expressible by Eq. 1, is systematically explored. We believe this type of model-based statistical approach to studying specific physical or biological models holds great promise.
The ability of KuraNet to generate computational hypotheses about collective motion in particularly complicated dynamical systems (Sec. Results III A) and to learn data-driven cluster synchrony patterns (Sec. Results III B), represents only the early steps of this modelbased approach. For example, note that our pairwise parallel processing of nodes relies on the fact that node features are statistically independent. Future work could relax this assumption and introduce more sophisticated coupling functions capable of learning more complicated relations among nodes in line with current work on complex network generation [39]. Similarly, future image segmentation work could use more hierarchical network structures like those in the standard convolutional neural networks [40] and in earlier oscillatory models [36]. Hierarchical networks whose couplings represent invariant features of natural images are likely necessary for training an oscillatory model to synchronize multiple images across a data set (as opposed to one at a time). Furthermore, our method could be applied to dynami-cal regimes beyond synchrony, for example the elusive glassiness speculated to appear in some oscillator networks [23,24]. Importantly, the framework we have called "KuraNet" can be easily applied beyond its namesake if Eq. 1 is replaced with another dynamics of interest. In that light, the proposed method becomes a way of casting the more general problem of understanding collective motion under conditions of quenched disorder in the framework of statistical learning. For example, the outlined technique could be used to learn coupling functions in neural systems or epidemiological networks, where disorder takes the form of heterogeneous cell biophysics [41] and infection properties [42], respectively. These applications will have to adapt the basic KuraNet approach to domain idiosyncracies, and we are intrigued by how this adaptation may give rise to insights in both the natural sciences and machine learning.