Higher-order interactions improve optimal collective dynamics on networks

Collective behavior plays a key role in the function of a wide range of physical, biological, and neurological systems where empirical evidence has recently uncovered the prevalence of higher-order interactions, i.e., structures that represent interactions between more than just two individual units, in complex network structures. Here, we study the optimization of collective behavior in networks with higher-order interactions encoded in clique complexes. Our approach involves adapting the Synchrony Alignment Function framework to a new composite Laplacian matrix that encodes multi-order interactions including, e.g., both dyadic and triadic couplings. We show that as higher-order coupling interactions are equitably strengthened, so that overall coupling is conserved, the optimal collective behavior improves. We find that this phenomenon stems from the broadening of a composite Laplacian's eigenvalue spectrum, which improves the optimal collective behavior and widens the range of possible behaviors. Moreover, we find in constrained optimization scenarios that a nontrivial, ideal balance between the relative strengths of pair-wise and higher-order interactions leads to the strongest collective behavior supported by a network. This work provides insight into how systems balance interactions of different types to optimize or broaden their dynamical range of behavior, especially for self-regulating systems like the brain.


I. INTRODUCTION
Complex networks provide the structural architecture for dynamical processes from a wide array of disciplines, and therefore their study constitutes an important fundamental area of research in physics, mathematics, biology, and engineering [1][2][3].Collective behaviors, i.e., consensus and synchronization, play particularly critical roles in the functionality of systems in many applications, with recent interest paid to applications including brain oscillations [4][5][6], cell signaling [7,8], and power grids [9,10].Moreover, various combinations of local dynamics with different microscopic and macroscopic topological network properties have been shown to give rise to a wide range of novel collective behaviors, including explosive synchronization transitions [11,12], chimera states [13,14], and macroscopic chaos [15,16], thus having important effects on system functions.
In addition to typical pairwise/dyadic interactions in network-coupled systems, recent work points to the presence of higher-order, e.g., triadic, interactions in both brain networks [17][18][19][20][21] and generic limit-cycle oscillator systems [22,23].The presence of such interactions is often encoded in simplicial complexes or hypergraphs [24][25][26][27] and has prompted the network science community to develop tools to better understand the impact of such higher-order interactions on collective dynamics.To date, a handful of studies have explored the role of higher-order interactions in collective dynamics in heterogeneous systems [28][29][30][31][32][33][34][35], but unlike real naturallyoccurring or engineered systems that are often optimized for a particular task, these initial studies tend to utilize random configurations or mean-field assumptions.At present, the role of higher-order interactions for optimized systems is largely * persebastian.skardal@trincoll.eduunknown for collective dynamics and other dynamical processes.
In this paper, we study collective dynamics in networks with higher-order interactions, focusing on collective behavior in optimized systems.To quantify the optimal collective behavior supported by a given network structure with higher-order interactions, we introduce a composite Laplacian matrix, which encodes the collective dynamics and network structure at multiple orders in a weighted simplicial complex and generalizes the Synchrony Alignment Function (SAF) framework [36] to this case.For the case of simple dyadic interactions, the SAF has been used to uncover the critical properties needed to optimize collective behavior in networks with heterogeneous dynamics and has proven to be flexibly adaptable to a wide range of realistic constraints and scenarios [37][38][39][40][41], as it encodes the interplay between heterogeneous dynamical units and heterogeneous network structure.We emphasize that in this context optimal refers to a system being as strongly synchronized as possible.
Generalizing the SAF framework, applying it in this new context, and analyzing the spectral properties of the composite Laplacian reveals important new properties of networks with higher-order interactions.Specifically, as higher-order interactions are strengthened in a system at the expense of weakening dyadic interactions to conserve the total coupling in the network, the eigenvalue spectrum of a composite Laplacian broadens.This, in turn, increases the dominant eigenvalue, which is key to improving the optimal state supported by the network.Complementing this improvement of the optimal collective state, the broadening of the eigenvalue spectrum also increases the overall range of possible states.This phenomenon contrasts sharply with synchronization of identical oscillators, where optimization of a network for identical synchronization reduces to a contraction of the eigenvalue spectrum [42,43].We close by exploring a realistic constrained optimization problem where local dynamics are not freely tunable, but must be allocated from a pre-defined set, and it is revealed that a network's ideal configuration is realized by a nontrivial, critical balance between the strength of dyadic and triadic interactions.

II. DYNAMICS AND MODELING
We begin with a higher-order generalization of the Kuramoto model [32,44] that consists of N phase oscillators whose states θ i , for i = 1, . . ., N , evolve according to Here, ω i is the natural frequency of oscillator i, K 1 and K 2 are coupling strengths that are associated with 1-and 2-simplex interactions, respectively, A is a 1-simplex adjacency matrix, and B is a 2-simplex adjacency tensor.We assume the network to be unweighted and undirected so that A ij = A ji = 1 if and only if a link exists between oscillators i and j, and and only if a triadic interaction exists between oscillators i, j and l.
While the 1-and 2-simplex coupling topologies may in general be uncorrelated for the case of a general hypergraph, here we assume the system corresponds to a simplicial complex so that the existence of a triadic interaction (i, j, l) requires the existence of dyadic interactions (i, j), (j, l) and (l, i).(Formally, the boundary of any 2-simplex in the simplicial complex must also be contained in the simplicial complex.)Moreover, we restrict our attention here to 'clique complexes' [45] in which all triangles give rise to 2-simplices, which allows the 3-tensor B to be completely determined by matrix A, i.e., The respective coupling strengths in Eq. ( 1) are scaled by the 1-and 2-simplex mean degrees k (1) and k (2) , which are population averages of the 1-and 2-simplex degrees k This scaling ensures that the overall connectivity is maintained between the 1-and 2-simplex structure.In other words, by conserving the sum K = K 1 + K 2 , we fix the overall amount of coupling in the network, regardless of the specific topologies encoded in A and B. To this end, we introduce a new bias parameter α ∈ [0, 1], defined via K 1 = (1 − α)K and K 2 = αK, so that α ≈ 0 corresponds to a system where 1simplex interactions are stronger than 2-simplex interactions and vice-versa if α ≈ 1.
Examples of a small toy network with 1-and 2-simplex dominated coupling are illustrated in Fig. 1 (a) and (b), respectively, where dyadic and triadic interactions are shaded to denote relative interaction strengths.Lastly, we note that other higher-order interaction terms may exist in other formulations of a higher-order Kuramoto model [22,23].We find that these yield qualitatively similar results as what is presented below, and so we focus our attention on the combination of 1-and 2simplex interactions in Eq. ( 1), and address additional triadic coupling terms in Appendix A.

III. OPTIMIZED SYSTEMS
Since our focus here is on optimization, which again we emphasize refers to a system being as strongly synchronized as possible, we consider the strongly synchronized regime where |θ j − θ i | ≪ 1, allowing us to linearize Eq. (1) to or in vector form, where L = (1 − α)L (1) + αL (2) is a composite Laplacian that is a weighted average of the first and second-order Laplacians, which we define, respectively as L (1) = (D (1) − A (1) )/ k (1)  and (2) .The matrix L (1) is simply a scaled version of the typical combinatorial Laplacian with D (1) = diag(k N ) and A (1) = A, while L (2) encodes the 2-simplex interactions with N ) and A (2) = A * (A 2 ) T , where * represents the Hadamard (i.e., element-wise) product.In addition to serving as the linear approximation of the nonlinear dynamics given in Eq. ( 1), Eqs. ( 2) and (3) also describe a forced consensus dynamics on the same network structure with higherorder interactions.Optimizing the synchronization dynamics using the linear approximation is equivalent to optimizing the consensus dynamics.
To optimize Eqs. ( 1)-( 3) we enter the rotating reference frame θ → θ + ω t (which allows us to effectively set the mean frequency to zero in both the nonlinear and linear dynamics), and we search for fixed points.Applying the Moore-Penrose pseudoinverse of the composite Laplacian [46], are the eigenvalues of L and its eigenvectors {v j } N j=1 form an orthonormal basis for R N , yields the fixed point From the viewpoint of consensus dynamics, the degree of consensus may be evaluated directly by the variance of the fixed point, θ * 2 /N .On the other hand, the degree of synchronization in the higher-order Kuramoto model is given by the magnitude r of the order parameter z = re iψ = N −1 N j=1 e iθj , which represents the centroid of all oscillators when placed on the complex unit circle.To leading order, the degree of synchronization of the fixed point is r ≈ 1 − θ * 2 /2N .Thus, consensus and synchronization dynamics are both optimized by minimizing the variance of the fixed point, θ * 2 /N .Using the form of L † given above and that θ * 2 = θ * , θ * , we have that The function J(ω, L) is known as the Synchrony Alignment Function (SAF), which was first introduced in Ref. [36] in the context of an objective function for optimizing the synchronization properties of a network of heterogeneous oscillators.Minimizing J(ω) serves to optimize θ * 2 /N and r and can be explored under a wide variety of constraints [37][38][39][40][41]. Inspecting the contributions to the SAF, we note that each term corresponds to a squared projection of the frequency vector ω onto the eigenvector v j that is scaled by inverse square of the associated eigenvalue λ j .Thus, under the constraint of fixing the variance the frequency vector to σ 2 , the collective behavior is strengthened by aligning the frequency vector ω as close as possible with the most dominant eigenvectors (those associated with larger eigenvalues) and orthogonalizing ω as best as possible to the least dominant eigenvectors (those associated with smaller eigenvalues).Thus, the optimal solution is obtained by setting ω = σ √ N v N .

IV. HIGHER-ORDER INTERACTIONS IMPROVE COLLECTIVE BEHAVIOR FOR OPTIMIZED SYSTEMS
In our first experiment, we highlight that random and optimized systems generally behave very differently, especially in the context of higher-order interactions.Specifically, we will show for optimized systems that collective behavior is improved by a stronger reliance on higher-order interactions, whereas it is diminished for random systems.Since simplicial complexes are geometrically embedded [47] we consider a class of noisy geometric networks [48] that contain both geometrically constrained and geometrically unconstrained edges between nodes uniformly placed on the unit disc in R 2 .With connected triangles, i.e., 2-simplexes, arising from geometrically constrained edges, we tune the prevalence of triadic interactions using a probability p ∈ [0, 1]: (i) with probability p each of the total M = N k (1) /2 edges is placed between the two closest nodes that are not yet connected and (ii) with probability (1 − p) each edge is placed randomly, where k (1) is the target mean 1-simplex degree.Thus, p tunes the prevalence of low-dimensional geometry in the network: in the limit p → 1 the network is purely geometric, while in the limit p → 0 the network is Erdős-Rényi [49].In Appendix B we provide a more complete algorithm implementing the network model described above.
Taking one such network of size N = 500 with mean degree k (1) = 10 and p = 0.25, we illustrate the effect of higher-order interactions on optimizing collective dynamics in Fig. 2(a).We plot the synchronization error 1 − r from direct simulations of Eq. ( 1) as a function of K for four cases, all under the constraint that the natural frequency vector has unit variance.First we consider the fully 1-simplex dominated case, i.e., α = 0 so that coupling is purely dyadic, and plot the results for random and optimal choices of natural frequencies in open and closed red triangles, respectively.Note that the optimal choice of natural frequencies outperforms the random case, given by a set of natural frequencies drawn from the standard normal distribution, by about an order of magnitude.Next, we set α = 0.8, thereby strengthening higher-order interactions at the expense of pair-wise interactions, and plot the results for random and optimal choices of natural frequencies in open and closed blue circles, respectively.We note here that all simulations are done using Heun's method with a time step of ∆t = 0.02, intregrating over a transient of 5 × 10 3 timesteps and then averaged over a steady state of 2×10 3 time steps.We also plot the predicted synchronization error, given by J(ω, L)/2K 2 , for each case in dashed curves, which accurately capture the dynamics for sufficiently large coupling.
This example highlights a critical feature of higher-order interactions in networks and their effect on collective dynamics.In particular, focusing on the optimal cases, the presence of higher-order interactions improves the optimal collective behavior supported by the system.Moreover, this phenomenon is generic: the more 2-simplex dominated a network is (i.e., the larger α is), the better the optimal states become.This is illustrated in Fig. 2(b), where over an ensemble of 10 3 networks built using the same parameters as the network using in Fig. 2(a) we plot the value of the SAF as a function of the bias parameter α for randomly chosen frequencies and the optimal choice in open and closed squares, respectively.(The average over this ensemble is plotted with dashed curves indicating one standard deviation up and down.)Specifically, we see that as α increases, thus making the the network more 2-simplex dominated, the optimal state improves very smoothly and monotonically while the random states worsen.Thus, strengthening higher-order interactions in collective network dynamics not only improves the optimal states, but also widens the range of possible states that are supported.

V. BROADENING OF COMPOSITE LAPLACIAN EIGENSPECTRUM UNDERLIES DICHOTOMY FOR OPTIMIZED AND NON-OPTIMIZED SYSTEMS
To explain and further illustrate the improvement that occurs in collective network dynamics as a result of increased higher-order interactions, we investigate the spectral properties of the composite Laplacian L = (1 − α)L (1) + αL (2) .Importantly, from Eq. ( 5) we can see that while the structure of the eigenvectors of L dictate the geometry of the optimal choice for the frequency vector ω, it is the eigenvalues that give insight into the quality of these optimal states.As an example, in Fig. 3(a) we plot the eigenvalue spectrum of L averaged across 10 3 networks of size N = 500 and built using the model described above with mean degree k (1) = 10 and p = 0.25 for α = 0 (solid blue), 0.4 (dashed red), and 0.8 (dot-dashed green).Note that as α increases and the higherorder interactions strengthen at the expense of pairwise interactions, the eigenvalue spectrum becomes broader.
In fact, it is the broadening of the eigenvalue spectrum, and specifically the increase in the dominant eigenvalue λ N , that corresponds to improving the optimal states, since, given the optimal choice ω = σ √ N v N , we have θ * 2 /N = J(ω, L)/K 2 = σ 2 /(Kλ N ) 2 .Here, we provide rigorous analytical insight on this mechanism by computing exactly the mean and variance of the eigenvalue spectrum in terms of moments of the various degrees using the trace of different powers of L. First, due to the conservation of the overall weighting of L (1) and L (2) , the mean is always conserved to one: i / k (2) ] = 1.Next, the variance Var(λ) = λ 2 − λ 2 = N −1 Tr(L 2 ) − N −2 Tr 2 (L) of the eigenvalue spectrum about this mean is given by where [See Appendix C for the derivation of Eq. ( 6).]In particular, varying α interpolates the variance between k (1)2 / k (1) 2 + 1/ k (1) − 1 and k (2)2 / k (2) 2 + q /(4 k (2) 2 ) − 1 in the extremes where connections are completely dominated by 1-simplex and 2-simplex coupling, respectively.Thus, when the latter form of the variance is larger, which we may expect when the 2-simplex degree distribution is more heterogeneous than the traditional 1-simplex degree distribution, strengthening higher-order interactions in turn broadens the eigenvalue spectrum of L. In general, as α is varied the eigenvalues interpolate between their respective values for L (1) and L (2) , however their intermediate behavior is more complicated and left for future research.
In Fig. 3(b), we plot the mean variance of the spectral density as a function of α which we calculated from the same ensemble as in panel (a) (indicating standard deviation with dashed curves).We observe a monotonic increase in the variance of the eigenvalue spectrum as higher-order interactions are strengthened, which is consistent with the broadening shown in panel (a).The extremal eigenvalues λ 2 and λ N follow this trend, decreasing and increasing, respectively, as illustrated in Figs.3(c) and (d).Moreover, we show over a full range of networks, from completely random to strongly geometric, that the 2-simplex degree distribution does in fact tend to be more heterogeneous than the 1-simplex degree distribution, thereby yielding improved collective dynamics as higher-order interactions are strengthened.In panel (e), we plot the 2-vs 1-simplex degrees for a single realization of the networks described above for parameters p = 0.01, 0.5, and 1, representing random, partially geometric, and completely geometric cases.The concave-up trend for each case suggests that the 2-simplex degree distribution is in fact more heterogeneous than the 1-simplex degree distribution.For a more concrete picture, we plot in panel (f) the quantities k (1)2 / k (1) 2 and k (2)2 / k (2) 2 (in blue circles and red crosses, respectively) for the network model discussed above across a full range of the parameter p, representing completely random networks (p ≈ 0) to completely geometric networks (p ≈ 1).Each data point represent the mean over an ensemble of 10 3 networks, with dashed curves representing one standard deviation.Here we see explicitly that over the full range we have that generically k (2)2 / k (2) 2 > k (1)2 / k (1) 2 , indicating that the phenomenon by which higher-order interactions improve optimal collective network dynamics in fact holds over a broad family of both random and geometric networks.
Furthermore, we may use this spectral analysis to shed light on more than just the optimal states, but also the worst possible state and random cases more broadly.First, analogous to the manner in which making the system more 2-simplex dominated increases λ N , and in turn promotes optimal collective dynamics, the complementary decrease in λ 2 results in poorer worst-case collective dynamics, which would result in setting the frequency vector proportional to the first non-trivial eigenvector, ω ∝ v 2 .Moreover, random frequency arrangements may be understood as follows.Constraining ω to unit variance σ 2 , it may be expanded using the eigenvector basis of L, ω = N j=2 c j v j , with N j=2 c 2 j = N σ 2 .Since heterogeneities are random and independent of network structure, the expected value of each coefficient is E[c j ] = ± N/(N − 1) thus and the expected value of J(ω, L) is given by where the average is taken over all eigenvalues except for the trivial eigenvalue λ 1 = 0.When higher-order interactions are then strengthened, broadening the eigenvalue spectrum, the decrease of the smaller eigenvalues λ 2 , λ 3 , . . ., which tend towards zero [See Fig. 3(a)] has a stronger effect on λ −2 than the increase of the larger eigenvalues . . ., λ N −1 , λ N .The overall effect of broadening the eigenvalue distribution by strengthening higher-order interactions is then increasing the expected value of the SAF and poorer expected collective behavior, even though, as we have seen above, the optimal states are improved.

VI. GEOMETRIC CONSISTENCY OF OPTIMAL SOLUTIONS
Next we explore the robustness of optimal solutions as the bias parameter is varied.In particular, we find a geometric consistency in the optimal and near-optimal choices of ω across a range of α.To illustrate this phenomenon, we consider the optimal choices of ω * for a handful of given bias parameter values α * = 0, 0.2, . . ., 1 and plot in Fig. 4(a) these values of the SAF J(ω * , L) as a function of α, with choices α * = 0 and 1 plotted in blue and yellow, respectively (and intermediate values interpolating these colors), for an ensemble of 10 3 networks of size N = 100 with mean degree k (1) = 10 and p = 0.25.
First, we note that, although one may expect the minimum of these curves to occur at α = α * (denoted by the open circles), the improvement of optimal collective behavior as α increases causes this minimum to occur at another value α > α * (denoted by the closed circles).Thus, after optimizing a system at a particular bias α = α * , increasing α results in improved collective behavior even without redesigning the frequency vector.Second, notice that for the case α * = 1 (given by the yellow curve) the value of the SAF J(ω * , L) remains quite close to the true optimal J(ω, L) for all other α.We see that this also generalizes: the SAF J(ω * , L) for a given value of α * remains close to the optimal SAF J(ω, L) for all α < α * .Thus, an optimal solution found for a 2-simplex dominated networks remains consistently near-optimal as the network becomes more 1-simplex dominated.
This phenomenon is due to a particular geometric property whereby the near-optimal subspaces of R N for two different values of α, namely the subspaces spanned by the eigenvec- tors v j of L with large eigenvalues λ j for the two values of α are largely overlapped.This is illustrated in Fig. 4(b) for α = 0.2 and α = 0.8, where we plot the logarithm (base-10) of the squared projection v j (0.2), v i (0.8) 2 for each (i, j) pair, using the same ensemble of networks as above and denoting v j (α) as the j th eigenvector of L for a bias parameter α.Note the strong diagonal feature indicating a strong alignment of v j (0.2) with v i (0.8) for the corresponding eigenvector i = j and other nearly adjacent eigenvectors i ≈ j.This geometric consistency likely has to do with structure of clique complexes, namely that the 2-simplex structure is defined precisely by the 1-simplex structure and all non-zero entries of the matrix A (2) come from non-zero entries of A (1) , and likely do not exist for more general cases of hypergraphs and simplicial complexes that are not clique complexes where 1-and 2-simplex structures may be uncorrelated.

VII. CONSTRAINED OPTIMIZATION
Lastly, we consider optimization of collective behavior in a more constrained scenario.Rather than just constraining the variance of a frequency vector ω and allowing frequencies to be freely chosen otherwise (thereby allowing them to be aligned perfectly with a particular eigenvector), we assume a randomly chosen initial frequency vector is given and may only be modified by a perturbation of constrained size.Denoting this perturbation by δω = ω new − ω, we then constrain the relative size δω / ω while maintaining the variance of the frequency vector itself.This perturbation may be optimally designed in terms of the eigenvector expansion ω = N j=2 c j v j by orthogonalizing away from the eigenvectors with smallest associated eigenvalues in order to eliminate the largest contributions to the SAF.To do this, we let δω = N j=2 β j v j and, for as large k as possible, let β j = −c j for j = 2, . . ., k, β j = c j for j = k + 1, . . ., N − 1, and , resulting in ω new = N j=2 γ j v j with γ j = 0 for j = 2, . . ., k, γ j = c j for j = k+1, . . ., N −2, and γ N = c N 1 + k j=2 c 2 j /c 2 N .Note that this both orthogonalizes ω new against the eigenvectors with smallest eigenvalues while increasing the alignment with v N in order to conserve the variance of ω.In Fig. 5(a), we plot the resulting SAF J(ω, L) averaged over 10 2 networks from an ensemble using the same parameters as in Figs. 2 and 3 after imposing such a perturbation of sizes δω / ω = 0, 0.4, 0.8, and 1.2 (plotted in blue circles, red triangles, green crosses, and black squares) to a random frequency vector with normally distributed entries.Note that the maximum possible perturbation that conserves the standard deviation of the frequencies is δω / ω = 2, obtained by ω new = −ω.
Another realistic possibility is that, before a perturbation is applied, the frequencies are (nearly) optimally rearranged to obtain a permutation of the initially given frequency vector.Here we obtain such a permutation using a simple acceptreject algorithm that interchanges randomly chosen pairs of frequencies if the exchange decreases the SAF.This preprocessing technique allows for more efficient perturbations, as we see in Fig. 5(b), which improve upon the results in Fig. 5(a).In particular, in such a constrained optimization scenario, we observe that there is often an ideal balance of dyadic to triadic interactions, i.e., a critical value of α that lies between zero and one, for a given perturbation size.
This phenomenon can be viewed as a combination between the cases of random frequencies, where higher-order interactions impede collective, and the case of optimal (freelytunable) frequencies, where higher-order interactions improve collective behavior.Specifically, the presence of a constraint allows higher-order interactions to improve the constrained optimal states, but only to a certain point since the purely optimal choice of frequencies is unattainable, as the frequency vector cannot be precisely aligned with the eigenvector vector v N .

VIII. DISCUSSION
Given the role of collective behavior in the function of physical, biological, and neurological systems where higherorder interactions may play a critical role in shaping system dynamics, understanding how higher-order interactions balance with dyadic interactions and affect collective behavior is an important question for a wide range of disciplines and applications.In this paper, we have addressed the topic of optimization for collective behavior in networks with higherorder interactions, focusing on clique complexes, and found that as higher-order interactions are equitably strengthened relative to dyadic interactions, optimal collective behavior improves.This phenomenon stems from the broadening of the eigenvalue spectrum of a composite Laplacian matrix that encodes the collective dynamics and network structure at multiple orders and generalizes the Synchrony Alignment Function framework to this important case.In particular, as the spectrum broadens the dominant eigenvalue(s) increase, which leads to this improvement.Moreover, we find that optimal solutions are robust over different balances between the relative strengths of dyadic and triadic interactions and that the broadening of the eigenvalue spectrum also widens the range of possible collective states supported by the network.We also find in more tightly constrained optimization scenarios that an ideal balanced between dyadic and triadic interactions occurs at a nontrivial, critical value of the bias parameter for networks to support the strongest possible collective behavior.Interestingly, the improvement of optimal collective behavior stemming from the broadening of the eigenvalue spectrum for the case of heterogeneous dynamical units lies in contrast to the case of identical units, where optimal networks stem from the concentration of the non-trivial eigenvalue spectrum [42].
These results shed light on the question of why higher-order interactions may be important in various applications that exhibit collective behavior.In particular, by modifying one's balance between dyadic and higher-order interactions, a system may self-regulate not only to modify (improve or worsen) its current collective state, but also broaden or contract the set of possible collective states that the system may support under further modification of the individual unit's local dynamics (here given by the oscillators' natural frequencies).This observation may be particularly useful for generating hypotheses for systems that exhibit a strong tendency for re-organization and self-regulation such as the brain where empirical evidence suggests that higher-order interactions play a role in collective behavior [17][18][19][20][21] and an optimal range of dynamic behavior is crucial for function [50,51].
This work introduced a composite Laplacian matrix that encodes network structure at multiple orders to generalize the SAF framework for optimizing collective behavior [36].The SAF framework has already been generalized for a number of scenarios involving optimization, including optimization with directed interactions [37], finding optimal network perturbations [38], optimizing networks with chaotic oscillators [39], addressing uncertainty in local dynamics [40], and uncovering geometric unfolding of networks [41].Future work may address many of these generalizations in the context of networks with higher-order interactions.Moreover, since we have restricted our attention to the case of clique complexes, it remains an open question how our findings would change if one were to consider other network structures such as more general simplicial complexes, hypergraphs, and non-geometric networks.

FIG. 1 .
FIG.1.Weighted simplicial complexes encode the balancing of multi-order interactions.An illustration of a small network with (a) 1-vs (b) 2-simplex dominated coupling.Shading indicates the relative strength of dyadic and triadic interactions which after rescaling by the respective mean degree, conserve the total coupling strength such that the bias parameter α equitably tunes 1-vs 2-simplex interactions.

FIG. 2 .
FIG.2.Optimal synchronization in networks with higher-order interactions.(a) The synchronization error 1 − r vs K for random (open symbols) and optimal (closed symbols) frequencies for two choices of the bias parameter: α = 0 (red triangles) and 0.8 (blue circles), representing cases where interactions are exclusively defined by 1simplexes and dominated by 2-simplexes, respectively, for a noisy geometric network (see text).(b) The Synchrony Alignment Function (SAF) J(ω, L) as a function of α for randomly-allocated (open squares) and optimal (closed squares) frequencies averages over 10 3 networks.

FIG. 3 .
FIG.3.Spectral properties of a composite Laplacian.(a) The eigenvalue spectrum P (λ) of the composite Laplacian L for α = 0 (solid blue), 0.4 (dashed red), and 0.8 (dot-dashed green) obtained from 10 3 networks of size N = 500 with mean degree k = 10 and p = 0.25.(b) The variance of the eigenvalue spectrum along with the extremal eigenvalues (c) λ2 and (d) λN from the same ensemble.(e) 2-simplex degrees k(2) vs 1-simplex degrees k(1) for a single network realization and (f) the quantities k (1)2 / k (1) 2 and k (2)2 / k (2) 2 (blue circles and red crosses, respectively) obtained from an ensemble of 10 3 networks as a function of the parameter p.

2 FIG. 4 .
FIG. 4. Geometric consistency in spectral properties of the composite Laplacian.(a) For networks optimized at α * = 0, 0.2, . . ., 1 (blue to yellow), the SAF J(ω * , L) of this solution as a function of α given the optimal frequency vector ω * .Points along the curves at α = α * and the local minimum are denoted in open and closed circles, respectively.(b) The logarithm (base-10) of the squared projections v j (0.2), v i (0.8) 2 .Results are obtained from an ensemble of 10 3 networks of size N = 100 with mean degree k (1) = 10 and p = 0.25.

FIG. 5 .
FIG.5.Constrained optimization.The SAF J(ω, L) as a function of α obtained after optimal perturbations of sizes δω / ω = 0 (blue circles), 0.4 (red triangles), 0.8 (green crosses), and 1.2 (black squares) are applied to a randomly drawn vector of frequencies (a) without and (b) with preprocessing the frequency vector using a (near) optimal permutation.Results are obtained from an ensemble of 10 2 networks of size N = 100 with mean degree k (1) = 10 and p = 0.25.