Network dynamics on graphops

Interacting particle systems, or more generally, dynamical systems on graphs/networks, form one of the major building blocks of modern science [1, 2]. Within the last three decades, they have permeated virtually all disciplines, ranging from molecular scales [3] to neuroscience [4], systems biology [5], machine learning [6], social science [7], epidemiology [8], and transportation networks [9] up to climate science scales [10]. For all-to-all coupled systems, lattice systems, and various special classes of network dynamics with similar dynamics at each vertex, there is a well-developed theory to pass to a mean-field approximation in many theoretical frameworks [11–16]. One considers the limit of an infinite graph [17]

$$\begin{equation*}n\to \infty \end{equation*}$$

so the particle/vertex number tends to infinity. If one is interested in the probability density ρ = ρ(u, t) to find a vertex in a certain dynamical state $u\in {\mathbb{R}}^{k}$ at time t, one may often derive a differential equation for ρ. One common form for k = 1 is

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho V\left(\rho \right)\right),\end{equation} \tag{ 1 }$$

where the map V can be derived from the dynamics of each vertex [11, 18]. If each vertex also is influenced by Gaussian noise, a second-order term ∂_uu(⋅) commonly appears as well in (1). The equation (1) is sometimes referred to as Liouville equation, continuity equation [19], and/or Vlasov equation, while the second-order equation is more known as Fokker–Planck and/or Kolmogorov equation. Here we shall refer to this class of equations as Vlasov–Fokker–Planck equation (VFPEs).

Now it is evident that one would like to derive a far more general form of (1), which also takes into account cases, where the adjacency matrix $A={A}^{\left(n\right)}\in {\mathbb{R}}^{n{\times}n}$ of the interaction between finitely many vertices is not just the the matrix of a full graph nor a highly symmetric structure such as a lattice, i.e., one wants to derive mean-field limits for complex networks [20, 21]. There have been many recent steps forward to achieve this goal. Several mathematical approaches have been successful providing rigorous proofs for VFPEs, where nonlocal integral terms appear to take into account the heterogeneous coupling structure [14, 22]. However, a general theory, which would allow us to upgrade easily from particular cases or standard all-to-all mean-field limit VFPEs, to modern complex network structures is still lacking.

In this note, we propose that a path to achieve an elegant extension of all classical approaches is to incorporate an analytical approach to graph limit theory. Graph limit theory [23], i.e., taking n → ∞, has made significant progress in recent years. New limit structures for dense graphs via graphons [24], as well as for sparse graphs [25–27], have recently appeared. The theory is relatively technical and convergence notions are often hard to understand as they have relied on combinatorial structures [23]. Recently a new approach to unify and extend graph limit theory was proposed by Backhausz and Szegedy [28]. Their idea relies on viewing the adjacency matrix A⁽ⁿ⁾ from an operator-theoretic perspective. In the simplest setting, one takes a vector v in ${\mathbb{R}}^{n}$ and consider the 2 × n matrix formed from v and v^⊤A⁽ⁿ⁾. Sampling the columns of this matrix uniformly at random generates a measure μ_v. Hence, one may view graphs also via their associated measures. Different graphs, even of different sizes, can then be more easily compared as we only view them through their action; see section 2 for more details and a brief review of the class of graphops (graph operators), which we shall rely on. Backhausz and Szegedy show that in a suitably chosen metric, large classes of graphs do have non-trivial subsequential limits

$$\begin{equation*}{A}^{\left(n\right)}\to {A}^{\left(\infty \right)}\quad \text{as} n\to \infty ,\end{equation*}$$

where A^(∞) can also be viewed as an operator, but now acting as an abstract operator on a function space. We are going to show (formally) that exactly this viewpoint is the missing ingredient to start a general theory of VFPEs, and related classes of evolution equations, as limits of finite-dimensional network dynamics.

We are going to illustrate our reasoning primarily on a benchmark class of interacting particle systems, the famous Kuramoto model of coupled oscillators [13, 19]. This model is widely studied and has served as a benchmark case to understand self-organizing and adaptive nonlinear synchronization. It can be embedded into a more general class of kinetic-type models. We shall see, how these kinetic models also fit elegantly together with graphops. This approach really seems to break a barrier that has held back the application of tools from dynamical systems, functional analysis, and evolution equations to broad classes of large-scale network dynamics problems.

The paper is structured as follows: in section 2, we review the relevant parts of graphop theory. In section 3, we show how this theory also provides an elegant framework to deal with network dynamics on a finite-dimensional level, i.e., for finite graphs. For example, we show that the graphop viewpoint provides a completely natural way to intertwine in one set of equations the dynamics on the network and the correlation structure between network elements. Then we proceed to formally develop, why graphops do generalize previous approaches to VFPEs in section 4. We conclude in section 5 with an outlook toward future opportunities and challenges.

We build upon the theory of graphops from [28]. Let $\left({\Omega},\mathcal{F},\mu \right)$ be a probability space and let ${L}^{p}\left({\Omega},\mathbb{R}\right)={L}^{p}\left({\Omega}\right)$ denote the usual Lebesgue space for p ∈ [1, ∞]. A linear operator A : L^∞(Ω) → L¹(Ω) is called P-operator if the operator norm

$$\begin{equation*}{\Vert}A{{\Vert}}_{\infty \to 1}{:=}{\mathrm{sup}}_{v\in {L}^{\infty }\left({\Omega}\right)}\frac{{\Vert}Av{{\Vert}}_{1}}{{\Vert}v{{\Vert}}_{\infty }}\end{equation*}$$

is finite. Key examples of P-operators are matrices, which we want to view as adjacency matrices of graphs. Indeed, consider Ω = {1, 2, ..., n} ≕ [n] for $n\in \mathbb{N}$ with the uniform measure μ_[n] and sigma-algebra given by the power set 2^[n]. Then a vector $v\in {\mathbb{R}}^{n}$ defines a function $v:\left[n\right]\to \mathbb{R}$ by v(j) = v_j for j ∈ [n]. On ([n], 2^[n], μ_[n]) any matrix $A\in {\mathbb{R}}^{n{\times}n}$ is a P-operator acting on (row) vectors via the usual rule

$$\begin{equation*}\left(vA\right)\left(j\right)=\sum _{k=1}^{n}v\left(k\right){A}_{kj}=\sum _{k=1}^{n}{v}_{k}{A}_{kj}.\end{equation*}$$

Obviously, we can identify L^∞([n]) and L¹([n]) with ${\mathbb{R}}^{n}$ and finite-dimensional linear operators given by matrices are bounded in the associated operator norm.

A key advantage of general P-operators are their convergence properties, when one considers sequences of these operators via profiles. As an example, we consider again $A\in {\mathbb{R}}^{n{\times}n}$ and any (row) vector $v\in {\mathbb{R}}^{n}$, then we have $vA\in {\mathbb{R}}^{1{\times}n}$ so that we may form a matrix $M\in {\mathbb{R}}^{2{\times}n}$ with rows v and vA. Now we sample columns of M uniformly, which yields a probability measure μ_v on ${\mathbb{R}}^{2}$. The (1-)profile of A is given by the collection of measures $\left\{{\mu }_{v}:v\in {\mathbb{R}}^{2}\right\}$. One may generalize the notion of a 1-profile to a k-profile ${\mathcal{S}}_{k}\left(A\right)$ for a general P-operator where one uses k vectors v₁, v₂, ..., v_k as the first k rows and v₁A, v₂A, ..., v_kA as k further rows. Hence, the matrix M becomes a matrix of size ${\mathbb{R}}^{2k{\times}n}$; wlog one may restrict v_j ∈ [−1, 1]ⁿ as the results do not change upon scaling vectors. A random column in M yields a probability measure on ${\mathbb{R}}^{2k}$ and ${\mathcal{S}}_{k}\left(A\right)$ is the collection of all such measures.

For a general P-operator A, consider functions v₁, ..., v_k, v₁A, ...v_kA and let

$$\begin{equation}\mathcal{D}\left({v}_{1},{v}_{2},\dots ,\;{v}_{k},{v}_{1}A,{v}_{2}A,\dots ,\;{v}_{k}A\right)\end{equation} \tag{ 2 }$$

be the joint distribution of the 2k-tuple, which can be viewed as the pushforward ${\left({T}_{k,A}\right)}_{{\ast}}\mu$ of the measure μ under the map

$$\begin{equation*}{T}_{k,A}\left(x\right)=\left({v}_{1}\left(x\right),{v}_{2}\left(x\right),\dots ,\;{v}_{k}\left(x\right),\left({v}_{1}A\right)\left(x\right),\left({v}_{2}A\right)\left(x\right),\dots ,\;\left({v}_{k}A\right)\left(x\right)\right),\end{equation*}$$

where x ∈ Ω. The k-profile ${\mathcal{S}}_{k}\left(A\right)$ is the set of all probability measures of the form (2), where we go through all possible k-tuples of functions in L^∞(Ω, [−1, 1]). In summary, the two crucial ideas are to view more analytically graphs as operators and to use profiles to compare matrices/graphs/P-operators regardless of their size. Let A, B be two P-operators. One defines a metric

$$\begin{equation}{d}_{M}\left(A,B\right){:=}\sum _{k=1}^{\infty }{2}^{-k}{d}_{\text{H}}\left({\mathcal{S}}_{k}\left(A\right),{\mathcal{S}}_{k}\left(B\right)\right),\end{equation} \tag{ 3 }$$

where d_H is the Hausdorff distance for $X,Y\in \mathcal{P}\left({R}^{k}\right)$ (i.e., X, Y are subsets of probability measures on ${\mathbb{R}}^{k}$)

$$\begin{equation}{d}_{\text{H}}\left(X,Y\right){:=}\mathrm{max}\left\{{\mathrm{sup}}_{x\in X}\,{\mathrm{inf}}_{y\in Y}\,{d}_{\text{LP}}\left(x,y\right),{\mathrm{sup}}_{y\in Y}\,{\mathrm{inf}}_{x\in X}\,{d}_{\text{LP}}\left(x,y\right)\right\}\end{equation} \tag{ 4 }$$

and d_LP is the Lévy–Prokhorov metric. This metric is actually quite natural and given for ${\nu }_{1},{\nu }_{2}\in {\mathbb{R}}^{k}$ by

$$\begin{equation}{d}_{\text{LP}}\left({\nu }_{1},{\nu }_{2}\right){:=}\mathrm{inf}\left\{\varepsilon { >}0:{\nu }_{1}\left(U\right){\leqslant}{\nu }_{2}\left({U}^{\varepsilon }\right)+\varepsilon \quad \;\text{and}\;\quad {\nu }_{2}\left(U\right){\leqslant}{\nu }_{1}\left({U}^{\varepsilon }\right)+\varepsilon \quad \forall U\in {\mathcal{B}}_{k}\right\},\end{equation} \tag{ 5 }$$

where ${\mathcal{B}}_{k}$ is the Borel sigma-algebra on ${\mathbb{R}}^{k}$ and U^ɛ is the set of points having distance less than ɛ from U. It is well-known that convergence in the Lévy–Prokhorov metric essentially is equivalent to weak convergence of measures so to actually ensure the existence of non-trivial limits, it is well-balanced choice between retaining some complexity in the limit, yet still having some form of mean-field object.

The definition of P-operators generalizes in a natural way to ||A||_p→q. One may metrize the space of P-operators via

$$\begin{equation}{d}_{M}\left(A,B\right){:=}\sum _{k=1}^{\infty }{2}^{-k}{d}_{\text{H}}\left({\mathcal{S}}_{k}\left(A\right),{\mathcal{S}}_{k}\left(B\right)\right),\end{equation} \tag{ 6 }$$

and call the associated convergence notion action convergence. Then one can prove [28] that any sequence ${\left\{{A}^{\left(j\right)}\right\}}_{j=1}^{\infty }$ with uniformly bounded norm ||⋅||_p→q for p ∈ [1, ∞) and q ∈ [1, ∞] has a limit [28, theorem 2.14]. This convergence is seen to generalize the theory of graphons and even graphings (classical objects used for very sparse graphs) as well as many intermediate cases. For example, if we have a graphon [23] $W:{\Omega}{\times}{\Omega}\to \mathbb{R}$ with finite norm

$$\begin{equation}{\Vert}W{{\Vert}}_{p}\,{:=}\,{\int }_{{\Omega}{\times}{\Omega}}W\left(x,y\right)\;\mathrm{d}{\mu }^{2},\quad p\in \left[1,\infty \right),\;q\,{:=}\,\frac{p}{p-1},\end{equation} \tag{ 7 }$$

then we can define an associated P-operator by

$$\begin{equation}\left(f{A}_{W}\right)\left(x\right)\,{:=}\,{\int }_{{\Omega}}W\left(y,x\right)f\left(y\right)\;\mathrm{d}\mu .\end{equation} \tag{ 8 }$$

It is easy to see that the norm ${{\Vert}{A}_{W}{\Vert}}_{p\to q}$ is finite in this case and A_W is just an operator-theoretic viewpoint on the L^p-graphon W. Again, compactness results hold for the case of ||⋅||_p→q so that we get the existence of a limiting graphon/P-operator.

There are also special classes of P-operators, which are particularly relevant. One important class are graphops, which are positivity preserving and self-adjoint P-operators, which behave effectively like undirected graphs. This makes sense since in the finite-dimensional case, adjacency matrices of undirected graphs are graphops.

Once we have graphops, it is natural to ask, how these operators are going to appear in dynamical systems on networks. It is the main theme of this paper to understand their appearance in dynamics on a formal level (rigorous proofs seem to be out of reach in full generality at this point but it is expected that eventually such proofs will be available).

We always denote the dynamical state of a vertex by u_k(t) for k ∈ [n] and denote the current time by t ∈ [0, T) with some fixed T > 0.

3.1. The Kuramoto model

The classical Kuramoto model considers oscillators on a circle so ${u}_{k}\left(t\right)\in {\mathbb{S}}^{1}{:=}\mathbb{R}/\left(0\sim 2\pi \right)$. The dynamics is given by

$$\begin{equation}\frac{\mathrm{d}{u}_{k}}{\mathrm{d}t}=:\,{u}_{k}^{\prime }={\omega }_{k}+\frac{{p}_{0}}{n}\sum _{j=0}^{n}\mathrm{sin}\left({u}_{j}-{u}_{k}\right),\end{equation} \tag{ 9 }$$

where ω_k are given internal frequencies of the individual oscillators and p₀ ⩾ 0 is a parameter controlling the coupling. The Kuramoto model on complex networks is usually written as follows [29]

$$\begin{equation}{u}_{k}^{\prime }={\omega }_{k}+p\sum _{j=0}^{n}{A}_{kj}\;\mathrm{sin}\left({u}_{j}-{u}_{k}\right),\end{equation} \tag{ 10 }$$

where p ⩾ 0 is a parameter controlling the coupling, which may contain a certain scaling in n depending upon the type of graph defined by the adjacency matrix $A={\left({A}_{kj}\right)}_{k,j\in \left[n\right]}$. Of course, on a finite-dimensional level the interpretation of adjacency matrices as graphops on ([n], 2^[n], μ_[n]), as discussed in section 2, applies. So if we consider a matrix of phase differences

$$\begin{equation*}{v}_{jk}{:=}\mathrm{sin}\left({u}_{j}-{u}_{k}\right),\quad V{:=}{\left({v}_{jk}\right)}_{j,k\in \left[n\right]},\end{equation*}$$

this yields upon inserting it into the Kuramoto model

$$\begin{equation}{u}_{k}^{\prime }={\omega }_{k}+p\;\mathrm{diag}\left(AV\right),\end{equation} \tag{ 11 }$$

where diag(M) is the vector obtained from the elements on the diagonal of the matrix M. Note that (V, AV) contains the same information as

$$\begin{equation*}\left({v}_{\cdot 1}^{\top },{v}_{\cdot 2}^{\top },\dots ,\;{v}_{\cdot n}^{\top },{v}_{\cdot 1}^{\top }A,{v}_{\cdot 2}^{\top }A,\dots ,\;{v}_{\cdot n}^{\top }A\right)\end{equation*}$$

where we used A = A^⊤ as our graphs are undirected. Hence, (V, AV) can be viewed as an n-profile if we uniformly sample from it. Since the matrix of phase differences V does evolve in time, we see that the input to the dynamics of each oscillator is (in addition to its intrinsic frequency) driven by moving through a subset of the n-profiles. Even beyond identifying the driving, we can go further by setting

$$\begin{equation*}\omega =\left({\omega }_{1},\dots ,\;{\omega }_{n}\right),\quad u=\left({u}_{1},\dots ,\;{u}_{n}\right),\quad d=\left({\left(AV\right)}_{11},\dots ,\;{\left(AV\right)}_{nn}\right),\end{equation*}$$

to re-write our differential equations as u' = ω + pd, so right-multiplication and time-differentiation of the 1-profile (u, uA) gives

$$\begin{equation*}\left(u,\left(uA\right)\right)=\left({u}^{\prime },{\left(uA\right)}^{\prime }\right)=\left(\omega +pd,\omega A,+pdA\right).\end{equation*}$$

Therefore, this defines an evolution equation for the pair (u, uA), which yields an evolution equation on the space of measures ${\left({T}_{1,A}\right)}_{{\ast}}\mu$ also written by

$$\begin{equation}\mathcal{D}\left(u,uA\right)=\frac{1}{n}\sum _{j=1}^{n}{\delta }_{{u}_{j},{\left(uA\right)}_{j}}\end{equation} \tag{ 12 }$$

which are associated to a 1-profile. In summary, the dynamics of the Kuramoto model on finite graphs induces an evolution equation on measures, which is somewhat different from using the classical empirical measure

$$\begin{equation*}\frac{1}{n}\sum _{j=1}^{n}{\delta }_{{u}_{j}},\end{equation*}$$

which keeps only track of the positions. Indeed, the 1-profile also contains the action of A on u given by uA and thereby the relevant correlation structure embodied by the joint empirical distribution (aka the 'JEDi') given by (12). The importance of this novel viewpoint of the dynamics is that once we are on a graph, then we have to intertwine structure and dynamics [30], thereby capturing how vertices respond to connectivity. We will see that this theme re-appears on other modeling scales and for VFPEs below.

3.2. The Cucker–Smale model

The Cucker–Smale model [31] is one well-known model for swarming. One variant of it reads

$$\begin{equation*}{u}_{k}^{\prime \prime }=\frac{1}{n}\sum _{j=1}^{n}\left({u}_{j}^{\prime }-{u}_{k}^{\prime }\right)\psi \left(\vert {u}_{k}-{u}_{j}\vert \right),\end{equation*}$$

where u_k = u_k(t) is the position of agent k at time t and ψ is a given function regulating the type of communication between different agents. More generally, one may pose the model on a network via

$$\begin{equation*}{u}_{k}^{\prime \prime }=\frac{1}{n}\sum _{j=1}^{n}{a}_{kj}\left({u}_{j}^{\prime }-{u}_{k}^{\prime }\right)\psi \left(\vert {u}_{k}-{u}_{j}\vert \right).\end{equation*}$$

We set ã_kj := a_kjψ(|u_k − u_j|), and obtain by re-writing the model in vectorized form as a first-order system and by right-multiplication with A, the following set of equations

$$\begin{equation}\begin{aligned}\hfill {u}^{\prime }\hfill & =v,\hfill \\ \hfill {v}^{\prime }\hfill & \hfill =\frac{1}{n}\left[\tilde{A}v-\text{Diag}\left(\tilde{A}{\tilde{A}}^{\top }\right)u\right],\hfill \\ \hfill {\left(uA\right)}^{\prime }\hfill & \hfill =vA,\hfill \\ \hfill {\left(vA\right)}^{\prime }\hfill & \hfill =\frac{1}{n}\left[\tilde{A}\left(vA\right)-\text{Diag}\left(\tilde{A}{\tilde{A}}^{\top }\right)\left(uA\right)\right],\hfill \end{aligned}\end{equation} \tag{ 13 }$$

where Diag(M) denotes the matrix obtained from M by only keeping the entries on the diagonal and setting all other entries to zero. Of course, the equations for (uA, vA) still depend directly on values of u via the definition of Ã. It is interesting to see that we effectively obtain now an evolution of measures via the 2-profile (u, v, uA, vA). The structure of the equations entails that the second-order nature of the model transfers to profiles, i.e., there is a constraint on the 2-profile in the same sense as for usual second-order ODEs.

3.3. Kinetic models

Instead of particular models, one can also look at broader classes of interacting systems frequently employed in kinetic theory [11, 16, 32]. A form commonly found in the literature is:

$$\begin{equation}{u}_{k}^{\prime }=\sum _{j\ne k}f\left({u}_{j},{u}_{k}\right).\end{equation} \tag{ 14 }$$

Evidently, one can extend this to a model with a complex network structure by writing

$$\begin{equation}{u}_{k}^{\prime }=\sum _{j=1}^{n}{a}_{kj}f\left({u}_{j},{u}_{k}\right).\end{equation} \tag{ 15 }$$

In this general abstract form, if we set f(u_j, u_k)  ≕  f_jk = f_jk(u) we get with a matrix F(u) = (f_jk(u)) that

$$\begin{equation*}{u}^{\prime }=\mathrm{diag}\left(AF\left(u\right)\right)\quad \text{and}\quad {\left(uA\right)}^{\prime }=\mathrm{diag}\left(AF\left(u\right)\right)A,\end{equation*}$$

as the evolution equation via the 1-profile (u, uA). Of course, this formulation is far too general to see any special structure regarding the evolution equation induced on the measures contained in profiles. If we linearize the evolution equations around a steady state u_* we obtain

$$\begin{equation*}{U}^{\prime }={D}_{u}{\left[\mathrm{diag}AF\left(u\right)\right]}_{u={u}_{{\ast}}}U\quad \text{and}\quad {\left(UA\right)}^{\prime }={D}_{u}{\left[\mathrm{diag}AF\left(u\right)\right]}_{u={u}_{{\ast}}}UA.\end{equation*}$$

So locally near steady states the evolution equation for the 1-profile are identical for both components, regardless of the precise underlying kinetic/particle model. This already hints at the fact, that profiles and the operator-theoretic viewpoint via graphops is well-suited for nonlinear dynamics and functional analysis techniques. This viewpoint will re-appear again for VFPEs below.

Having observed that graphops and profiles can provide a very elegant re-interpretation as well as augmentation of finite-dimensional network dynamics models, we now proceed to formally take limits n → ∞ to obtain 'mean-field' models. We consider the mesoscale, in the sense that we are no longer interested in the state of individual vertices on a graph but only in the probability distribution of the state of a vertex.

4.1. The Kuramoto model

We have already seen in section 3.1 that the classical Kuramoto model provides a good starting point. It is well-studied, and its homogeneous mean-field limiting equation is well-known. We recall its classical formal derivation [13, 19] here because this derivation will be important for us to reflect back upon in the context of graphops below.

4.1.1. All-to-all coupling

Starting from (9), suppose we have very large number of oscillators and their intrinsic frequencies are distributed according to a density g = g(ω) with g(ω) = g(−ω). One introduces the complex order parameter as

$$\begin{equation}r{\text{e}}^{\text{i}\psi }{:=}\frac{1}{n}\sum _{j=1}^{n}{\text{e}}^{\text{i}{u}_{j}}.\end{equation} \tag{ 16 }$$

Multiplying the order parameter by ${\text{e}}^{\text{i}{u}_{k}}$, and taking imaginary parts, one easily checks that now the Kuramoto model (9) can be re-written as

$$\begin{equation}{u}_{k}^{\prime }={\omega }_{k}+{p}_{0}r\;\mathrm{sin}\left(\psi -{u}_{k}\right).\end{equation} \tag{ 17 }$$

In particular, the kth oscillator feels all other oscillators via a single mean-field parameter. Next, let

$$\begin{equation*}\rho \left(u,t,\omega \right)\;\mathrm{d}u\end{equation*}$$

denote the fraction of oscillators with frequency ω between u and u + du at time t. Then by construction ρ is a probability density

$$\begin{equation*}{\int }_{0}^{2\pi }\rho \left(u,t,\omega \right)\;\mathrm{d}u=1,\quad \rho {\geqslant}0.\end{equation*}$$

In the limit n → ∞, we can formally re-write the order parameter (16) as

$$\begin{equation}r\left(t\right){\text{e}}^{\text{i}\psi \left(t\right)}={\int }_{0}^{2\pi }{\int }_{\mathbb{R}}{\text{e}}^{\text{i}u}\rho \left(u,t,\omega \right)g\left(\omega \right)\;\mathrm{d}\omega \;\mathrm{d}u.\end{equation} \tag{ 18 }$$

Indeed, the last formula can be derived as a limit of the sum in (16) using the law of large numbers [19]. Yet, it is very crucial to note that this law of large numbers argument only gives us a mean-field as it produces only the mean of the order parameter completely discarding correlations. Furthermore, the approach implicitly assumed that each random variable has finite variance. Yet, if these assumptions hold then we know that the resulting Liouville/continuity equation for the probability density ρ for the Kuramoto model should be given by

$$\begin{equation*}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho v\right),\quad v\left(u,t,\omega \right){:=}\omega +{p}_{0}r\;\mathrm{sin}\left(\psi -u\right).\end{equation*}$$

This equation can be re-written using the definition of r(t) in (18), multiplying by a suitable exponential and taking the imaginary part as

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho \left(\omega +{p}_{0}{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}\mathrm{sin}\left(\tilde{u}-u\right)\rho \left(\tilde{u},t,\tilde {\omega }\right)g\left(\tilde {\omega }\right)\mathrm{d}\tilde {\omega }\mathrm{d}\tilde{u}\right)\right).\end{equation} \tag{ 19 }$$

The previous argumentation can be made fully rigorous to derive the mean-field Vlasov-type equation (19). Yet, the argument is highly non-trivial and the currently most elegant way of proof proceeds via a so-called Dobrushin-type bound [33, 34], which is effectively a Gronwall estimate on a metric space of measures, usually carried out in a Wasserstein metric [11, 35, 36]. We remark that if one would consider the Kuramoto model with stochastic perturbations, then one is going to obtain a (nonlinear) Fokker–Planck equation [18], which appears in many classical stochastic coupled oscillator models, where the coupling appears through a mean field, such as the Desai–Zwanzig model [37]. Hence, one sometimes refers to (19) as a Vlasov–Fokker–Planck (VFPE) equation.

4.1.2. Complex network heuristics

Having dealt with the classical all-to-all coupled VFPE case, it is natural to think about extensions to treat the Kuramoto model (10) on complex networks. Recall that we have assumed that the underlying graph is undirected and unweighted so that the adjacency matrix $A={\left({A}_{ij}\right)}_{i,j\in \left[n\right]}$ is symmetric and binary A_ij ∈ {0, 1} for all i, j ∈ [n]. There exists a known formal approach [29] to try to 'save' the classical mean-field argument from section 4.1.1. One defines a local order parameter as

$$\begin{equation*}{r}_{k}{\text{e}}^{\text{i}{\psi }_{k}}{:=}\sum _{j=1}^{n}{A}_{kj}{\text{e}}^{\text{i}{u}_{j}}.\end{equation*}$$

If the graph is sufficiently well-connected and effectively still behaves like an all-to-all coupled system, one is then tempted to make the ansatz of a single global field

$$\begin{equation*}r{\text{e}}^{\text{i}\psi }=\frac{1}{{\kappa }_{k}}{r}_{k}{\text{e}}^{\text{i}{\psi }_{k}},\end{equation*}$$

where κ_k is the degree of vertex k, i.e., the global field is locally proportional to the local field weighted by the degree. With this assumption one obtains that (17) can now be written as

$$\begin{equation*}{u}_{k}^{\prime }={\omega }_{k}+{p}_{0}r{\kappa }_{k}\;\mathrm{sin}\left(\psi -{u}_{k}\right).\end{equation*}$$

In particular, these steps motivate that one might want to use the local degree of a vertex as the next approximation criterion to derive a mean field [38]. Hence, one option is to consider a probability density ρ(u, t, ω, κ) of oscillators having phase u, frequency ω, and degree κ at time t with the usual conditions

$$\begin{equation*}{\int }_{0}^{2\pi }\rho \left(u,t,\omega ,\kappa \right)\;\mathrm{d}u=1,\quad \rho {\geqslant}0.\end{equation*}$$

Evidently one cannot really do much with standard tools as just defining ρ still does not lead to a mean-field VFPE. If one assumes that the network is uncorrelated and has a degree distribution d(κ), then the probability that an edge has its end a vertex of phase u, degree κ and frequency ω at time t is

$$\begin{equation*}\frac{\kappa d\left(\kappa \right)}{\langle \kappa \rangle }g\left(\omega \right)\rho \left(u,t,\omega ,\kappa \right)\end{equation*}$$

where ⟨κ⟩ is the average degree of a vertex in the graph. This suggests that it might be helpful to define the order parameter now as

$$\begin{equation*}r{\text{e}}^{\text{i}\psi }{:=}{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}{\int }_{0}^{\infty }\frac{\kappa d\left(\kappa \right)}{\langle \kappa \rangle }g\left(\omega \right)\rho \left(u,t,\omega ,\kappa \right)\;\mathrm{d}\omega \;\mathrm{d}\kappa \;\mathrm{d}u.\end{equation*}$$

Now the same trick as previously, multiplying by an exponential and taking imaginary parts, yields VFPE equation

$$\begin{equation*}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho \left(\omega +p\kappa r\;\mathrm{sin}\left(\psi -u\right)\right)\right).\end{equation*}$$

The right-hand side of the VFPE could again be expressed now as some integral. Indeed, one may formally write the equation for the evolution of the average phase in the limit n → ∞ as

$$\begin{equation}{u}^{\prime }=\omega +\frac{p\kappa }{\langle \kappa \rangle }{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}{\int }_{0}^{\infty }g\left(\tilde {\omega }\right)d\left(\tilde {\kappa }\right)\tilde {\kappa }\rho \left(\tilde{u},t,\tilde {\omega },\tilde {\kappa }\right)\mathrm{sin}\left(\tilde{u}-u\right)\;\mathrm{d}\tilde {\kappa }\;\mathrm{d}\tilde{u}\;\mathrm{d}\tilde {\omega },\end{equation} \tag{ 20 }$$

which arises as a formal replacement of the sum in the Kuramoto model (10). This suggests a closed VFPE in the form

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho \left(\omega +\frac{p}{\langle \kappa \rangle }{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}{\int }_{0}^{\infty }g\left(\tilde {\omega }\right)d\left(\tilde {\kappa }\right)\tilde {\kappa }\rho \left(\tilde{u},t,\tilde {\omega },\tilde {\kappa }\right)\mathrm{sin}\left(\tilde{u}-u\right)\;\mathrm{d}\tilde {\kappa }\;\mathrm{d}\tilde{u}\;\mathrm{d}\tilde {\omega }\right)\right).\end{equation} \tag{ 21 }$$

The last considerations already show an emerging theme that we shall exploit later on: there is an effective integral operator acting on the density ρ, which depends upon the graph structure, in the mean-field VFPE equation on complex networks. In comparison to the classical all-to-all coupled VFPE (19), we have replaced in (21)

$$\begin{equation*}\rho \left(\tilde{u},t,\tilde {\omega }\right)\;\,\text{by}\,\;{\int }_{0}^{\infty }\frac{d\left(\tilde {\kappa }\right)\tilde {\kappa }}{\langle \kappa \rangle }\rho \left(\tilde{u},t,\tilde {\omega },\tilde {\kappa }\right)\;\mathrm{d}\tilde {\kappa }.\end{equation*}$$

Of course, this replacement entails implicit assumptions, most prominently that the degree distribution alone is enough to approximate the dynamics, e.g., higher-order correlations between degrees are discarded. So far, it has been difficult to validate such assumptions mathematically beyond simulations, which shows that a more general abstract framework may be needed.

4.1.3. Kuramoto on graphons

Regarding the previous considerations in section 4.1.2 on integral operators, it is now no surprise, how the VFPE equation on graphons should look [22]. Let $W:{\Omega}{\times}{\Omega}\to \mathbb{R}$ be a graphon as defined in (7) and recall that it defines an associated P-operator A_W via (8). For graphons we know that we may identify Ω = [0, 1] with the unit interval. So for x ∈ Ω one now obtains the VFPE

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho \left(\omega +p{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}{\int }_{{\Omega}}g\left(\tilde {\omega }\right)W\left(x,\tilde {x}\right)\rho \left(\tilde{u},t,\tilde {\omega },\tilde {x}\right)\mathrm{sin}\left(\tilde{u}-u\right)\;\mathrm{d}\tilde {x}\;\mathrm{d}\tilde{u}\;\mathrm{d}\tilde {\omega }\right)\right)\end{equation} \tag{ 22 }$$

for a density ρ = ρ(u, t, ω, x) of oscillators having phase u, frequency ω, and graphon position x at time t. The position x ∈ [0, 1], where we evaluate a graphon has the natural interpretation in the graph limit as the position/index of a vertex, while W(x, y) provides a limiting version of the adjacency matrix, i.e., how x is connected to y. In comparison to the classical all-to-all coupled VFPE (19), we have replaced in (22)

$$\begin{equation*}\rho \left(\tilde{u},t,\tilde {\omega }\right)\;\,\text{by}\,\;{\int }_{{\Omega}}W\left(x,\tilde {x}\right)\rho \left(\tilde{u},t,\tilde {\omega },\tilde {x}\right)\;\mathrm{d}\tilde {x}.\end{equation*}$$

It is important to observe that (22) still relies crucially on the fact that the graph limit can be encoded by a single scalar-valued integral kernel W(x, y). This covers only a relatively small set in the space of all graph limits as explained in [28]. Yet, the generalization to far more general classes beyond graphons is becoming physically evident if one aims to remove the integral operator restriction.

4.1.4. Kuramoto on graphops

Let A : L^p(Ω) → L^q(Ω) be a graphop as defined in section 2. Suppose it arises as a limit from a sequence of adjacency matrices

$$\begin{equation*}A={A}^{\left(\infty \right)}={\mathrm{lim}}_{n\to \infty }\;{A}^{\left(n\right)},\end{equation*}$$

where convergence is understood with respect to the metric d_M. From the Kuramoto model on complex networks (10) one now formally obtains the VFPE

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho \left(\omega +p{\int }_{0}^{2\pi }{\int }_{\mathbb{R}}g\left(\tilde {\omega }\right)\left({A}^{\left(\infty \right)}\rho \right)\left(\tilde{u},t,\tilde {\omega },x\right)\mathrm{sin}\left(\tilde{u}-u\right)\;\mathrm{d}\tilde {x}\;\mathrm{d}\tilde{u}\;\mathrm{d}\tilde {\omega }\right)\right)\end{equation} \tag{ 23 }$$

for a density ρ = ρ(u, t, ω, x) of oscillators having phase u, frequency ω, and graphop coordinate x at time t. The position x ∈ Ω is now more abstract, yet it still has a meaning in the sense that it should present the position/index of a vertex. Instead of an integral operator representation, for a P-operator, we have to think of the action of a graph, which is actually the primary insight that recently arose in the graph theory literature as discussed in section 2. In (23), we see, how this insight can be carried over to easily writing abstract VFPE mean-field type equations, even if the graph is neither all-to-all, nor uncorrelated, nor homogeneous, nor dense, nor a special sparse graph. In comparison to the classical all-to-all coupled VFPE (19), we have abstractly replaced in (23)

$$\begin{equation*}\rho \left(\tilde{u},t,\tilde {\omega }\right)\;\,\text{by}\,\;\left({A}^{\left(\infty \right)}\rho \right)\left(\tilde{u},t,\tilde {\omega },\tilde {x}\right).\end{equation*}$$

Note that this formulation is much cleaner, more general, and much easier to comprehend in comparison to other approaches.

4.2. Kinetic models on graphops

The approach in the previous section for the Kuramoto model is now relatively easy to formally generalize to more abstract kinetic models such as (15). One simply takes the normal VFPE for (14), which would read

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho V\left(\rho \right)\right),\end{equation} \tag{ 24 }$$

where V is completely computable from f [11]. Then one lifts the action of the finite-dimensional adjacency matrices A⁽ⁿ⁾ onto a limiting graphop A^(∞), and replaces the density in the vector field driven part V of the VFPE

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho V\left({A}^{\left(\infty \right)}\rho \right)\right),\quad \rho =\rho \left(u,t,x\right),\;x\in {\Omega}.\end{equation} \tag{ 25 }$$

All the operator-theoretic properties of the finite-dimensional graphs A⁽ⁿ⁾, which persist in the limit A^(∞), actually can now be used to analyze the VFPE. We know from [28] that graphops do carry a lot of information via the graph limit for large classes of graphs, so we may strongly expect that an analysis of (25) can now proceed along classical lines of dynamical systems theory. For example, steady states ρ_* will satisfy

$$\begin{equation*}0={\partial }_{u}\left({\rho }_{{\ast}}V\left({A}^{\left(\infty \right)}{\rho }_{{\ast}}\right)\right).\end{equation*}$$

Linearization is possible via the normal Fredholm derivative in many cases. Suppose that we have an evolution equation

$$\begin{equation}{\partial }_{t}\rho =-{\partial }_{u}\left(\rho V\left({A}^{\left(\infty \right)}\rho \right)\right)\quad \rho =\rho \left(t\right)\in \mathcal{X},\end{equation} \tag{ 26 }$$

in a suitable Banach space $\mathcal{X}$ where the P-operator part is such that it defines a well-defined evolution in $\mathcal{X}$. Then the linearization is just formally given by the chain rule and product rule

$$\begin{equation}{\partial }_{t}{\Xi}=-{\partial }_{u}\left({\Xi}V\left({A}^{\left(\infty \right)}{\rho }_{{\ast}}\right)+{\text{D}}_{\rho }V\left({A}^{\left(\infty \right)}{\rho }_{{\ast}}\right){A}^{\left(\infty \right)}{\Xi}\right)=:\text{\,\,}{\mathcal{L}}_{\infty }{\Xi}\end{equation} \tag{ 27 }$$

for ${\Xi}={\Xi}\left(t\right)\in \mathcal{X}$. The linear operator ${\mathcal{L}}_{\infty }$ can now be analyzed using classical tools from functional analysis such as spectral theory. As in section 3.3 on the finite-dimensional level, the operator ${\mathcal{L}}_{\infty }$ carries in an intertwined way the information about

• (S1)  
  the shape of the steady state via ρ_*,
• (S2)  
  the linearization of the coupling function f via D_ρV,
• (S3)  
  the graph spectral information via A^(∞).

Hence, the graphop viewpoint provides in a completely natural way the opportunity to lift spectral information from the finite-dimensional setting to the graph limit in the context of stability analysis.

4.3. Remarks on applications

Applying the strategy (S1)–(S3) for concrete graphs is beyond the current work as it requires more detailed investigations of the spectral theory of graphops, which is currently under development. However, we would like to outline possible abstract steps one might take to link topology and dynamics similar to the general idea of the master stability function formalism [39] following up on the last section. In the master stability function framework, one usually makes directly assumptions about the shape or symmetry of steady states. On the mean-field level for graphops, it is much more natural to make assumptions about A^(∞)ρ_*. For simplicity, assume A^(∞) : L^p → L^q maps ρ_* to a finite-dimensional sub-space of L^q, so that we can express the image A^(∞)ρ_* as a linear combination

$$\begin{equation*}{A}^{\left(\infty \right)}{\rho }_{{\ast}}=\sum _{l=1}^{L}{\rho }_{{\ast},l}{e}_{l}\end{equation*}$$

for a suitable Schauder basis e_l and computable coefficients ρ_∗,l. Then we have to study, when

$$\begin{equation*}\lambda \text{Id}\left(\cdot \right)+{\partial }_{u}\left(\cdot V\left(\sum _{l=1}^{L}rh{o}_{{\ast},l}{e}_{l}\right)+{\text{D}}_{\rho }V\left(\sum _{l=1}^{L}{\rho }_{{\ast},l}{e}_{l}\right){A}^{\left(\infty \right)}\left(\cdot \right)\right)\end{equation*}$$

has no bounded inverse to get the elements $\lambda \in \mathbb{C}$ of the spectrum. Suppose we may understand this problem by studying the two operators

$$\begin{equation}{\partial }_{u}\left(\cdot V\left(\sum _{l=1}^{L}{\rho }_{{\ast},l}{e}_{l}\right)\right)\quad {\partial }_{u}\left({\text{D}}_{\rho }V\left(\sum _{l=1}^{L}{\rho }_{{\ast},l}{e}_{l}\right){A}^{\left(\infty \right)}\left(\cdot \right)\right)\end{equation} \tag{ 28 }$$

separately; note that this is already a (potentially) strong assumption similar to the case, when two matrices have to commute to easily determine the spectrum of their sum. The first operator in (28) does no longer depend on A^∞ as this dependence is in the finitely many coefficients ρ_∗,l. The second operator in (28) can be viewed as a composition of three operators, which we may abbreviate ∂_u D_ρV A^(∞). If these operators commute then we have

$$\begin{equation*}\text{spec}\left({\partial }_{u}\;{\text{D}}_{\rho }V\;{A}^{\left(\infty \right)}\right)\subseteq \text{spec}\left({\partial }_{u}\right)\text{spec}\left({\text{D}}_{\rho }V\right)\text{spec}\left({A}^{\left(\infty \right)}\right).\end{equation*}$$

In the commuting case, we then get that the spectrum of the graphop (encoding the topology) multiplied by the spectrum of the linearization of the vector field (encoding the dynamics) provide a spectral bound, which can then be translated into a stability bound for the state ρ_* as we just have to check, when the spectral bound ensures that the spectrum lies in the left half of the complex plane, as usual.

Although we did make quite strong assumptions in the last calculations, it shows that the general approach to link topology to network dynamics is still possible. More precisely, the approach shows that it is a next natural step to establish spectral properties of graphops and to check when A^(∞)ρ_* is practically computable.

We have shown that the viewpoint of a more abstract operator-theoretic approach is extremely elegant to provide a general theory for network dynamics, when we do not have a simple coupling structure. This applies to the finite-dimensional microscopic case as well as the limiting case for VFPEs.

Although we have physically derived, by using analogies to previous approaches in the mean-field context for homogeneous/all-to-all coupling, a suitable VFPE formulation on a limiting graphop, several questions remain. From a kinetic theory perspective, we would like to prove that the limiting VFPE indeed approximates over a finite time scale, the underlying finite-dimensional network dynamics for large n. We conjecture that similar proofs¹ via a Dobrushin-type estimate can also work in the graphop context, when the structure of the graph limit theory is taken into account.

From the dynamics perspective, it would be very interesting to study spectral theory for operators involving graphops as one element of their definition in more detail as we have indicated in section 4.3. A natural first example, where such a theory could be used is again the Kuramoto model on a complex network.

Furthermore, from the physical perspective it is important to gain a better understanding of graph limit procedures in the sense of particle interactions. We believe that the finite-dimensional observations we provided hint at a correlation structure viewpoint being particularly important to provide the correct physical interpretation. From the viewpoint of direct applications in various sciences, it seems reasonable to assume that testing examples and analyzing finite-size effects are going to be of particular practical importance.

C K acknowledges partial support of the DFG via the SFB/TR 109 Discretization in Geometry and Dynamics, partial support by the TiPES (Tipping Points in the Earth System) project funded the European Union's Horizon 2020 research and innovation program under grant agreement No. 820970 and support by a Lichtenberg Professorship of the VolkswagenStiftung. C K also thanks the TUM International Graduate School of Science and Engineering IGSSE for support via the project SEND (Synchronization in Evolutionary Network Dynamics). C K also appreciated general discussions with Sebastian Throm and Tobias Böhle regarding limits of network dynamics. The comments of two anonymous referees have helped to improve the presentation of the results.