Coupled Dynamics on Hypergraphs: Master Stability of Steady States and Synchronization

In the study of dynamical systems on networks/graphs, a key theme is how the network topology influences stability for steady states or synchronized states. Ideally, one would like to derive conditions for stability or instability that instead of microscopic details of the individual nodes/vertices rather make the influence of the coupling structure visible. On graphs, the master stability function is an important such tool. Here we generalize the master stability approach to hypergraphs. A hypergraph coupling structure is important as it allows us to take into account arbitrary higher-order interactions between nodes. As for instance in the theory of coupled map lattices, we study Laplace type interaction structures in detail. Since the spectral theory of Laplacians on hypergraphs is richer than on graphs, we see the possibility of new dynamical phenomena. More generally, our arguments provide a blueprint for how to generalize dynamical structures and results from graphs to hypergraphs.


INTRODUCTION
Dynamical systems on networks are a fundamental part of the theory of complex systems [BBV08,PG16]. A common situation in network dynamics is that one would like to infer dynamical conclusions just from the underlying network structure. This has led to the introduction of the master stability function formalism [PC98], see also the exposition in [New10]. The idea is to assume sufficient symmetry and/or common dynamics for each individual node/vertex, which then makes it possible to re-write stability conditions for steady states, or even more complicated synchronized solutions, in terms of network data. Examples of network data in this context are spectra, e.g., of the graph Laplacian or the adjacency matrix [Chu97]. The master stability function approach has been successfully applied in many applications, particularly in the context of synchronization of oscillators [BP02,ADGPV06,NMLH03]; see also the surveys [ADGK + 08, DB14]. However, just considering binary interactions modelled by a network/graph is often insufficient in applications. One then needs generalizations of graphs. A first natural generalization are simplicial complexes [Hat02]. Simplicial complexes have appeared in several applications, e.g., in protein classification [CMW + 15], in percolation models for statistical physics [BKZ19], in computational neuroscience [GGB16], in modelling dynamics of social peer pressure [HK20], or in epidemiology [IPBL19,MGA20].
More generally, these results are examples that higherorder interactions [BAR16,GBMSA17,SA19] are relevant between nodes/vertices, where we note that the study of higher-order interactions has already quite a long history, particularly in ecology [Abr83,BC94]. While simplicial complexes form a very convenient mathematical structure, they are also somewhat rigid as not all possible higher-order interactions are allowed. This led to an interest to study more general hypergraphs, e.g., for cellular networks [KHT09], for opinion formation [LN13], for epidemic spreading [BKS16], or for social network analysis [ZL10]. For instance, consider collaboration relations among scientists (see for instance [PDJ19]). We may have scientists A, B, C that coauthor a paper, and there may also exist a paper written by A and B without C, as well as single author papers by A and C, but no others. This would be modelled by a hypergraph with vertices A, B, C and hyperedges {A}, {C}, {A, B}, {A, B, C}. Neither a graph nor a simplicial complex would be adequate to capture this structure. Therefore, in this paper, we study dynamics on hypergraphs. We shall generalize the general tool of master stability functions from graphs to hypergraphs. In particular, we derive general conditions for the linear stability of synchronized dynamics. We then turn to the important special class of Laplace type interactions where we can apply the recently developed spectral theory for hypergraph Laplacians [JM19]. At the end, we provide an outlook how our framework could be used as a blueprint to systematically generalize dynamical aspects of graphs to hypergraphs.

SETTING: STABILITY FOR SYSTEMS OF ODES
We briefly recall linear stability theory for systems of ordinary differential equations to fix the notation and the main ideas. Let us consider a set of units i = 1, . . . , N , called nodes or vertices, in the sequel, that are dynamically interacting with each other. This leads to a system of differential equations, where we assume that the state variables x i could be vector-valued, We may then also write (1) in matrix form A solution x * of (1) is called linearly stable, or simply stable for short in the sequel, if any solution ǫ of the linearization or in the more abstract version corresponding to (2) converges to 0 for t → ∞. Here, is the vector with components ∂Fi(x * 1 ,...,x * N ) ∂x α j , α = 1, . . . , m, and therefore, in (3) there is an implicit sum over α. Linear stability is simply a condition on the Lyapunov exponents of the tensor ∂Fi(x * 1 ,...,x * N ) ∂xj (note that this tensor will in general depend on time t, since we are not assuming that x * (t) is constant. The stability condition then can be expressed in terms of a Lyapunov exponent (see for instance [Arn98]), lim sup t→∞ 1 t log e tDF(x * (t)) < 1.
There are two special cases that are of particular interest.
1. The solution x * is constant in time, that is, steady or stationary. This means that for each i, x * i (t) = x * i (0) is independent of time t. Such a stationary state simply satisfies For such a solution, the stability condition is simply (5).
2. The solution x * (t) represents a synchronized state. This means that it is independent of the vertex i, that is, x * i (t) = x * j (t) for all i and j, and all t. To make such a solution feasible, we should also assume that F i is the same for all i. For the stability of synchronization, we only need to require that any non-synchronized solution of (3) converges to 0 for t → ∞.
In the sequel, we shall only consider the second case. The first case succumbs to a similar, but easier analysis.

INTERACTION ON NETWORKS
We now consider the situation where a vertex i does not interact indiscriminately with all other vertices but only maintains interactions with a subset of vertices; those vertices are called the neighbors of i, and one writes j ∼ i when j is such a neighbor of i. When one considers network interactions, these interactions are assumed to be pairwise only. That means that we may be able to write the dynamical system (1) in the form Here, f i is a self-interaction term of i, whereas g ij stands for the pairwise interaction between i and j. In order to make the interaction structure more explicit, one also considers systems of the form (see also [New10]) or (see also [PC98]) where the (vector-valued) dynamical functions f , g, h no longer depend on the vertices. The focus then is on the interaction matrix A = (a ij ) i,j=1,...N . The neighborhood structure can be included in that matrix by stipulating that a ij = 0 unless j ∼ i. We consider (9), as the analysis of (8) is similar. The resulting stability condition has been referred to in the literature as master stability condition. If one wishes to make synchronized dynamics possible, one usually assumes that does not depend on i. In that case, a synchronized solution x * of (9) would satisfy The linear stability equation (4) for (9) at a solution x * is ( where Id always denotes the identity operator of suitable size, which is simply the n-dimensional identity matrix in the context of (12). When we assume that the coupling matrix A can be diagonalized (for instance, if it is symmetric, i.e., a ij = a ji for all i, j), we let its eigenvalues be µ k , k = 1, . . . , N . Since 1 N is the identity matrix, we can decompose (12) into the corresponding modes ǫ k , that is, When we assume (10), one of the eigenvectors of A is constant. Therefore, at a synchronized state x * , we obtain a mode ǫ 1 (t) with ǫ 1 i (t) = ǫ 1 j (t) for all i, j. The evolution of the mode therefore leaves the synchronization manifold invariant. Synchronization is stable when all other modes decay. Let us consider the case where f = h. Then (13) becomes The stability condition then is (see [JJ02]) that is, where as in (5), is the maximal Lyapunov exponent of f (at the particular solution x * , but in order to have a general criterion, we may take the supremum over all solutions). The inequality (16) now separates and relates the condition for the dynamical update f and the network connectivity nas encoded in the coupling matrix A and its eigenvalues. In the interesting case, we have ℓ f > 1, that is, the dynamics generated by f is unstable. But if the eigenvalues µ 2 , . . . , µ N lie between −2 and 0 and satisfy (16), synchronization may still be a stable state. Similar to [JJ02], we now consider the case where Here, 0 ≤ σ ≤ 1 is a parameter and is the normalized Laplace operator of the network (see for instance [Chu97,Jos14] for the theory, but note that the conventions employed here are somewhat different from those in these references). The eigenvalues of ∆ satisfy where the eigenfunction for λ 1 = 0 is constant. The stability condition (16) then becomes that is, by (20), Thus, we need at the same time a lower bound for the first nonzero eigenvalue and an upper bound for the largest eigenvalue. λ 2 is controlled from below by the so-called Cheeger inequality [AM85,Dod84] which quantifies the cohesion of the graph. λ 2 is largest when the graph is complete, and of course, a complete graph is more conducive to synchronized dynamics than a less coherent one. In particular, λ 2 = 0 precisely if the graph is disconnected, and for such a graph, we obviously cannot expect dynamics to synchronize. In fact, when the graph has more than one component, the dynamics could be synchronized on each component, but not necessarily between components. Let us consider the case of two components Γ 1 , Γ 2 . An eigenfunction for λ 2 = 0 then is constant on each component (with the weighted sum of the constants being zero). When ℓ f > 1, but (22) is satisfied now for λ 3 , then what we may call the generalized synchronization manifold, that is, the family of dynamical states that are synchronized inside the two components only, is stable against perturbations by other eigenstates. Analogously, of course, for more than two components.λ N = 2 holds precisely if the graph is bipartite, and in fact the gap 2 − λ N quantifies the deviation from bipartiteness [BJ13]. On a bipartite graph, antiphase oscillations are possible, and thus, there again is an obstacle to synchronization carried by the mode associated with λ N . That is why we need the upper bound. Given ℓ f and the topology of the underlying graph, (22) then tells us whether we can find a range of coupling strengths σ for which synchronized dynamics are stable.

INTERACTION ON HYPERGRAPHS
So far, we have essentially summarized or reformulated known results. In particular, in the preceding section, we have considered dynamics on a network where the dynamics at each vertex is coupled with the dynamics of its neighbors. The network thus corresponds to a graph with edges defined by the neighborhood relations. Thus, all relations are binary. When we also want to include higher order interactions, as in many empirical systems, we need an underlying structure that is more general than that of a graph. We need a hypergraph. A hypergraph has a set V of vertices i = 1, . . . , N and a set H ⊂ 2 V of hyperedges h = 1, . . . , M . Thus, each hyperdedge is a set of vertices h = {i h1 , . . . , i hm h } where m h is the number of vertices contained in the hyperedge h. We can then consider types of dynamics analogous to those in equations (7). These can be written as We note that the number of arguments of an interaction function g ih now depends on the size m h of the hyperedge h. When we linearize (23), we therefore need the N × M incidence matrix I := (I ih ) defined by We observe that, for each i and j, Therefore, (24) Returning to the general system (23), its linearized version at a solution x * then is For the stability of x * , we need to check as before whether ǫ(t) → 0 as t → ∞ for any solution of (25). After this general result, we now want to discuss the possibility and the stability of synchronized dynamics on hypergraphs. Again, we require (10). When we want to consider the analogue of (8) or (9) and again assume uniform interaction functions, these functions will now still on the size of the hyperedgs, as the number of their arguments varies with the size m of the underlying hyperedge. Thus, we have functions g m . When we have an interaction matrix A = a ih , the dynamics then are of the form When, for instance a ih = I ih , (25) becomes We thus see (24) in action. For instance, we may consider the case where g m (y 1 , . . . , y m ) is a normalized symmetric function of its entries, for instance Importantly, we can again consider a Laplacian type coupling. The corresponding hypergraph Laplacian was constructed in [JM19], and we shall recall some of its properties. First of all, we need a little bit more structure than that of a simple hypergraph. We need what is called a chemical hypergraph in [JM19], that is a collection of vertices i = 1, . . . , N and a collection of oriented hyperedges h = 1, . . . , M . An oriented hyperedge is a non-empty ordered subset (V h , W h ) of 2 V × 2 V . The vertices in V h and W h are called the inputs and outputs of h. Changing the orientation of h simply means replacing As before, the stability condition couples the Lyapunov exponent of the dynamical nonlinearity f , the structure of the hypergraph as encoded by the eigenvaluesλ k of∆, and the coupling parameter σ. Indeed, if we replace in (18) the usual graph Laplacian by the hypergraph Lapla-cian∆, then we get a stability condition |1 − σλ k |ℓ f < 1 for k = 1, . . . , N, Note carefully, that although we have 0 ≤λ 1 ≤λ 2 ≤ · · · ≤λ N , we do not have the same strong bounds as for the usual graph Laplacian as presented in (20). Yet, we can still re-write (29) as, whereλ min is the smallest non-zero eigenvalue. Even for a connected hypergraph,λ 2 need not be greater than 0. This, in fact, leads to an interesting class of dynamics. Let us assume thatλ 1 = · · · =λ k = 0, butλ k+1 Then the class of dynamics that belong to eigenstates of the Laplacian for the eigenvalueλ = 0 is stable. This class can be larger than the locally synchronized dynamics. For instance, consider a graph with three vertices 1, 2, 3 and a single hyperedge with V = {1}, W = {2, 3}. One eigenstate forλ = 0 is constant, but another one is given by u(1) = 1, u(2) = u(3) = 1 2 . This would correspond to a dynamical state x * with g(x * 2 ) = g(x * 3 ) = 1 2 g(x * 1 ), which would be stable under our conditions. That is, the dynamical activity at 1 is equally split into the activities at 2 and 3, as prescribed by the topology of the hypergraph. -Conversely, it may also happen that all eigenvalues of a hypergraph are positive. Take, for instance, again three vertices, and for each i a hyperedge h i with V hi = {i}, W hi = {i + 1, i + 2}, counting the vertices mod 3. Then all eigenvalues are positive, see [JM19], precluding the possibility of synchronized dynamics. Furthermore, another difference with the graph case is that 2 does not give an upper bound toλ N . In fact,λ N is equal to N in some cases and it is not known yet whether this is the largest possible value for λ N . Nevertheless, the geometrical meaning of the largest eigenvalue does not change. It is in fact known that, given a hypergraph Γ with largest eigenvalueλ N , theñ whereλ ′ N is the largest eigenvalue of a bipartite hypergraph that has the same number of hyperedges as Γ and also the same number of inputs and the same number of outputs in each hyperedge (catalysts are not included). Also, the equality holds if and only if Γ is bipartite.
In summary, we find that once the hypergraph Laplacian appears in the dynamics directly, one can still derive a master stability condition. But one has to be careful, e.g., in treating the dimension of the synchronization manifold as well as possible degenerate additional neutral modes associated to zero eigenvalues, which may appear on a linear level for the hypergraph Laplacian. In addition, it is clear that hypergraph coupling can shift the stability regions. This lends some interest to results for a particular model in the special case of simplicial complexes [IPBL19]. However, note that our master stability conditions only operate on the level of the linearization. The case of higher-order interactions and bifurcations, where nonlinearities matter even locally, is far more involved [KB20].
Finally, we point out that while in this article, we have considered time-continuous dynamics, our scheme also applies to time-discrete dynamics. For instance, one can study the phenomenon of the synchronization of chaos [Kan84] on analogues of coupled map lattices on hypergraphs.

CONCLUSION AND OUTLOOK
In this work we have shown how to extend the master stability function framework from graphs to hypergraphs. In particular, we noticed how the spectral properties of the hypergraph Laplacian enter the stability condition, and how this changes the statements we may make regarding the interplay between network topology and dynamics. For example, it is now possible that the upper bound on the largest eigenvalue grows significantly, while already the smallest eigenvalue can be bigger than zero. Conversely, even for connected hypergraphs, the multiplicity of the eigenvalue 0 can be larger than 1, and this leads to interesting new classes of dynamics that are more general than synchronization, but may still be locally stable under appropriate conditions. Furthermore, we found that the incidence matrix plays an important role in hypergraph dynamics, and it interacts in a non-trivial way with the master stability condition(s). We point out that the approach we have taken here provides a general strategy for lifting results about dynamics on graphs to hypergraphs. The key is to identify the steps where the adjacency matrix or the graph Laplacian play key roles, and then replace them with analogous hypergraph objects. The spectral theory of hypergraphs is richer than that of graphs, and that lead us to identify new classes of dynamics that are more general than synchronization but for which we can still derive stability conditions analogous to those for synchronized dynamics on graphs.