An algebraic approach to the Kuramoto model

We study the Kuramoto model with attractive sine coupling. We introduce a complex-valued matrix formulation whose argument coincides with the original Kuramoto dynamics. We derive an exact solution for the complex-valued model, which permits analytical insight into individual realizations of the Kuramoto model. The existence of a complex-valued form of the Kuramoto model provides a key demonstration that, in some cases, re-formulations of nonlinear dynamics in higher-order number fields may provide tractable analytical approaches.

The dynamics of networks with many nodes and connections poses difficulties for mathematical treatment. While linear dynamics on a network of interacting units is straightforward to handle with matrix algebra, once nonlinearity is introduced at individual units, the dynamics of the system becomes difficult to study analytically. Here, we consider the Kuramoto model (KM), a paradigmatic example of nonlinear network dynamics describing the synchronization of coupled oscillators [1]. This model provides a canonical description of synchronization in nature, where populations of interacting units (from neurons to Josephson junctions and fireflies [2]) coordinate the timing of their behavior in the absence of a central coordinator [3]. Because of its wide applicability throughout many systems, the KM is a central consideration in the study of nonlinear dynamics.
Analytical approaches to the KM started with the original work of Kuramoto, who introduced the model in [4]. By changing to a rotating coordinate frame and passing to a continuum description, Kuramoto provided a description of the population dynamics in the infinite limit and the transition to synchrony at a critical coupling strength. Since this pivotal work, much research has gone on to analyze the dynamics of KM. In recent years, several studies have focused on introducing complex versions of the KM [5,6], with the goal of providing some analytical insight; however, finding a complex KM that permits an exact solution for the evolution of this system has proven evasive.
In this Letter, we provide a complex-valued matrix formulation of the KM whose argument corresponds to the original KM. We derive an explicit analytical solution for the dynamics using the eigenspectrum of the adjacency matrix. We then use this solution to study synchronization in individual realizations of the model. The existence of a complex-valued version of the KM permitting an analytical solution provides a key demonstration for potential analytical approaches to complex nonlinear dynamics at finite scales. The complex-valued KM and its exact solution provide an example that, in some cases, reformulations of nonlinear dynamics in different number fields may provide opportunities for analytical descriptions of nonlinear systems.
We start with the standard Kuramoto model (KM) on a general network of N nodes: where θ i (t) ∈ [−π, π] is the state of oscillator i ∈ [1, N ] at time t, ω i is the intrinsic angular frequency, κ scales the coupling strength, and element a ij ∈ {0, 1} represents the connection between oscillators i and j. The standard sine coupling in the interaction term causes phases of two connected oscillators i and j to attract, to an extent depending on the homogeneous coupling strength κ. Here, we first consider the KM defined on an undirected ring graph G RG , where nodes are arranged on a one-dimensional ring with periodic boundary conditions and connected to their k neighbors in each direction (later considering Erdős-Rényi and Watts-Strogatz graphs). We also consider first the case of homogeneous intrinsic frequencies ω i = ω ∀ i ∈ [1, N ] (with inhomogeneous case in Supplement [7]). For k = N/2 , G RG corresponds to the fully connected graph on N nodes, K N . It is important to note that the approach described here is general to the specific graph model we consider. Starting with Equation (1), we can change to a rotating coordinate frame [8] and set ω = 0 without loss of generality. We then subtract an additional imaginary component in the interaction term: We note this expression now implies θ i ∈ C and requires a scaling of the coupling strength (γ = 2κ/π); as we will show, the resulting complex-valued system admits a solution whose argument agrees with the original KM.
Multiplying both sides of Equation (2) by i, we have Applying Euler's formula, we can rewrite the system as: This equation now results in the following matrix form: where θ ∈ R N is the state vector across nodes, and A ∈ {0, 1} N ×N is the adjacency matrix encoding the connections between nodes in the network. We can then utilize the fact that the inverse of the matrix diag[e −iθ ] is diag[e iθ ] to arrive at Now, using the fact that mentioned in [6], we have that and letting x = e iθ , we havė whose general solution is We can now utilize a diagonalization of the adjacency matrix, A = V DV −1 , to rewrite the solution as and the problem of computing the matrix exponential is now reduced to computing the standard exponential function on the diagonal entries of D. We note that the adjacency matrix A is always diagonalizable for the undirected graphs considered here. Finally, as noted above, Equation (10) can be directly related to the original KM. Let θ = θ re + iθ im be the decomposition of θ into the real and imaginary parts. Then we have θ re is thus the argument of the analytical solution x. In particular, we can take θ re ∈ [−π, π]. We will use this solution to compare with the numerical integration of Equation (1). What remains in deriving an expression for the resulting temporal dynamics is to study the eigenspectrum of A. To do this, we note the adjacency matrix of the ring graph G RG is by definition a circulant matrix [9], and we can use the circulant diagonalization theorem (CDT) to obtain its eigenspectrum analytically. The CDT states that all circulants c ij = circ(c j ), where circ(c j ) is the circulant matrix constructed from the generating vector c j , are diagonalized by the same unitary matrix U with components r, s ∈ [1, N ] and that the N eigenvalues are given by Using these expressions, we can then evaluate Equation (10) in terms of the eigenspectrum of A to obtain a fully analytical evaluation of x(t), which we can then compare to numerical integration of the original KM.
We can now compare the argument of the analytical expression Equation (10) with the result obtained by numerical integration of Equation (1). As a first example, we considered a small network of N = 3 nodes and k = 1, resulting in the complete graph K 3 . Figure 1 (top left) shows the considered network along with the state variables θ i at each node. Initial conditions θ i (0) were selected randomly with uniform distribution U(−π, π). Starting from these random initial conditions and a homogeneous coupling strength κ = 1, this small network synchronizes over the course of the one second simulation (cf. dotted blue vertical lines across bottom three panels, Fig. 1). Equations were integrated using a forward Euler method with fine temporal precision (dt = 0.001 s) and were compared with results from a Runge-Kutta method [10]. Additional simulations were run at very high temporal precision (up to dt = 10 −6 s) to ensure the numerical validity of these results. Simulation code to reproduce all figures in this work will be made availabe at our GitHub site (http://mullerlab.github.io).
We next compare the analytical and numerical results for larger networks of N = 200 nodes on a ring graph with k = 100. To create a synchronization transition during a one second simulation interval, we scaled the coupling strengths accordingly (κ = 6/N ). Figure 2 shows the results of this comparison between the numerical (top) and analytical (bottom) evaluations. The numerical and analytical evaluations exhibit similar dynamics, from the macroscopic synchronization to the specific trajectories of individual oscillators. Interestingly, a small fraction of the nodes in the numerical simulation remain counterphase to the rest of the population, while this behavior does not occur in the analytical version.
To understand this point further, we systematically studied synchronization in the KM. To quantify the extent of synchronization, we use a standard approach to study the sum over oscillators: where r(t) ∈ C, |r(t)| ∈ [0, 1] is the order parameter at time t, and θ j (t) ∈ [1, N ] is the phase of oscillator j at time t (n.b. θ j,re (t) in the case of the analytical expression Equation 10). Figure 3 shows the time-average order parameter over 1-second simulations. As κ increases in the KM, the order parameter begins at a low value (desynchronized state) and increases until approaching unity (synchronized state) at values of κ ranging from 1 to 10. As observed in Figure 2, the numerical and analytical versions of the KM exhibit very similar macroscopic synchronization dynamics. Lastly, we considered the numerical and analytical solutions of the KM on undirected Erdős-Rényi and Watts-Strogatz random graphs. For each realization of a random graph, we obtain a numerical estimate of the eigenspectrum of its adjacency matrix A to use in the analytical expression Equation (10). We first considered the KM on an Erdős-Rényi random graph (G ER ), which displays synchronization dynamics similar to those previously observed (Fig. 4, top). We then considered the KM on a Watts-Strogatz network (G W S ), which is defined as a ring graph G RG where each node is first connected to its k neighbors in each direction and each edge is rewired to another node with uniform probability q [11]. The KM on G W S displays non-trivial spatiotemporal dynamics before converging to the synchronized state (Fig. 4, bottom middle). Importantly, these spatiotemporal dynamics were also well described by the analytical expression introduced here (Fig. 4, bottom right).
In this work, we have introduced a complex-valued formulation of the KM whose argument corresponds to the original Kuramoto dynamics. This formulation permits an exact analytical solution for individual realizations of the KM. Here, we have first considered the case of homogeneous intrinsic frequency; however, this approach generalizes to the inhomogeneous case (Supplement [7]). We then compared the analytical version to numerical integration of the KM equations and found the analytical version displays similar dynamics, from the macroscopic synchronization behavior to the specific trajectories of individual oscillators.
Previous research, including an inspiring technical report written by van Mieghem [6], has studied complexvalued formulations of the KM [5,12]; however, no exact analytical expression was previously obtained. In particular, the expression derived in [6] was noted to hold only for the repulsive cosine variant of the KM, and the linear reformulation in [5] requires tuning a parameter to create the correspondence to the original KM. Thus, the results reported in this Letter, while representing only an initial study utilizing the analytical expression in Equation (10), to the best of our knowledge represent the first analytical version of the Kuramoto model.
Importantly, we emphasize that the analytical version introduced here is valid at finite scales and for individual realizations of the KM. This analytical version allows future mathematical study of synchronization dynamics in networks with many nodes and connections, potentially using new tools from spectral graph theory, in addition to allowing one to obtain the future state of the system at an arbitrary future moment without numerical integration of the differential equations defined on the network. Because the KM has been extensively studied both as a model for neural dynamics [8,13,14] and as a fundamental model for neural computation [15], these results open up several possibilities for studying the connections between network structure, nonlinear dynamics, and computation. Recurrent connections have previously been shown to produce powerful computations through nonlinear interactions [16]. The approach introduced here may have applications in understanding such recurrent interactions, which have been increasingly acknowledged to play an important and unexplained role in visual processing in the brain [17]. Understanding more clearly the connection between networks and computation thus may have implications for fields such as neuroscience and beyond.
In this work, we have specifically studied the KM defined on a ring graph G RG , whose highly regular structure permits analytical study of the eigenspectrum of its adjacency matrix. The regular structure of this graph means that its adjacency matrix belongs to the class of circulant matrices, whose eigenspectrum can be calculated analytically. Further, when k is maximal, such that the number of neighbors to which a node is connected equals the rest of the graph, G RG corresponds to K N , the complete graph on N nodes. The KM defined on K N , in turn, corresponds to the case of all-to-all connectivity first con-sidered by Kuramoto [4]. For these reasons, we chose to focus on the KM defined on G RG in this work. In previous work [18], however, we have introduced an operatorbased approach to the structure of random graphs. In future work, we aim to extend the present results to understand the connection between nonlinear dynamics and random graphs with various structural features.
The existence of a complex-valued formulation of the Kuramoto model that permits an analytical solution raises an important example in nonlinear dynamics. While Equation (1), which includes the sine coupling interaction term and the network adjacency matrix, appears analytically intractable, in this case a reformulation in the complex domain provides an algebraic approach to the Kuramoto dynamics which serve as a canonical description for synchronization phenomena in nature. While this re-formulation may not generalize beyond the Kuramoto model, the ability of this approach to provide insight in this case suggests that representations in number fields beyond R may represent opportunities for future insight.