Reconstructing higher-order interactions in coupled dynamical systems

Higher-order interactions play a key role for the operation and function of a complex system. However, how to identify them is still an open problem. Here, we propose a method to fully reconstruct the structural connectivity of a system of coupled dynamical units, identifying both pairwise and higher-order interactions from the system time evolution. Our method works for any dynamics, and allows the reconstruction of both hypergraphs and simplicial complexes, either undirected or directed, unweighted or weighted. With two concrete applications, we show how the method can help understanding the complexity of bacterial systems, or the microscopic mechanisms of interaction underlying coupled chaotic oscillators.


I. INTRODUCTION
Higher-order interactions are present in ecosystems, where the way two species interact can be influenced by a third species [1], in social systems, where interactions in groups of three or more individuals naturally occur [2], in the brain cortex [3], and in many other complex systems [4].Recent studies based on mathematical tools such as simplicial complexes [5,6] and hypergraphs [7] have already demonstrated that the dynamics in presence of higher-order interactions can be significantly different from that of systems where interactions are exclusively pairwise [2,[8][9][10][11].How to infer and model higher-order interactions is then crucial for understanding the dynamics and functioning of complex systems [12,13].While in complex networks, the reconstruction problem, also known as the inverse problem, i.e. determining the network from the dynamics of a system, has been dealt with different techniques [14], the question on how to infer connectivity in presence of higher-order interactions is still open.
Concerning reconstruction in complex networks, two different types of approaches, which target either the functional or the structural connectivity of the system, have been developed.In the first type of approaches, functional networks are typically constructed from the network temporal evolution by evaluating correlation, Granger causality or transfer entropy among the signals of the different network units, or using Bayesian inference methods [15][16][17].In the second type of approaches, the underlying structural connectivity of a network is obtained from the network response to external perturbations [18], from its synchronization with a copy with adaptive links [19,20], or from the solution of optimization problems based on measurements of node time series, when the functional form of the node dynamics is known [21][22][23].
The reconstruction problem in the presence of higher-order interactions is more convoluted.Recently, the fundamental distinction between higher-order mechanisms, i.e. the presence of higher-order terms in the microscopic structure of the interactions, and higher-order behaviors, i.e. the emergence of higher-order correlations in the dynamical behaviour of a system, has been pointed out [24].The relationship between these two is not trivial as higher-order behaviors do not necessarily rely on higher-order mechanisms.However, for their identification, techniques that go beyond pairwise statistics are required.For instance, information-theoretic approaches to study multivariate time series (of node activities) based on hypergraphs [25], higher-order predictability measures (such as generalizations of Granger causality and partial information decomposition) [26], or simplicial filtration procedures [27] have been proposed to extract important information on higher-order behaviours that otherwise would not be visible to standard, i.e., network-based, analysis tools.Higher-order behaviors, which are likely due to the presence of higher-order mechanisms, can be identified by recently introduced techniques for statistical validation able to detect overexpressed hyperlinks [28,29].Other statistical approaches to the problem are based on Bayesian methods, and have been used to construct hypergraphs directly from pairwise measurements (link activities), even in cases where the higher-order interactions are not explicitly encoded [30,31].Statistical inference and expectation maximization are also at the basis of a method recently developed to reconstruct higher-order mechanisms of interaction in simplicial SIS spreading and Ising Hamiltonians with two-and three-spin interactions [32].However, this method can only be applied to binary time series data produced by discrete two-state dynamical models.
In this work, we propose an optimization-based approach to infer the high-order structural connectivity of a complex system from its time evolution, which works in the case of the most general continuous-state dynamics, i.e. when node variables are not restricted to take binary values.Namely, we consider a system made by many dynamical units (nodes) coupled through pairwise and higher-order interactions; we assume that the local dynamics and the functional form of interactions are known [33] or identifiable [34,35], and we propose a method to extract the topology of such interactions by solving an optimization problem based on the measurement of the time evolution of the node variables.With two concrete applications, we show that the method can effectively reconstruct which nodes are interacting in pairs and which in groups of three or more nodes.

II. RECONSTRUCTION OF THE INTERACTIONS
As a general model of a dynamical system of N nodes coupled through pairwise and higher-order interactions, we consider the following set of equations: with i = 1, . . ., N .Here x i (t) ∈ R m is the state vector of unit i, f i : R m → R m is the nonlinear function describing the local dynamics at node i, while g (d) : R m×(d+1) → R m are the nonlinear functions of order d, modeling interactions in groups of d + 1 nodes, with d = 1, . . ., D. The topology of (d + 1)-body interactions is encoded in the tensor A (d) , while the parameters σ 1 , . . ., σ D > 0 tune the strengths of the coupling at each order.
Here we want to infer the complete structural connectivity of a dynamical system, which means we want to reconstruct, not only the entries of the adjacency matrix A (1) = {a (1) ij } from the knowledge of the evolution of the state variables x 1 (t), . . ., x N (t), but also the higher-order interactions encoded by the tensors In doing this, we assume that the functions f i and g (1) , g (2) . . ., g (D) , are known.This is a reasonable assumption as the local dynamics of many real-world complex systems, as well as the functional forms of their interactions, have been well identified.For instance, well-established mathematical models to describe the dynamics of neurons and synapses, or the growth of a biological species when in isolation, or when in interactions with other species, are available.In absence of such models, we assume instead that, prior to the structural connectivity reconstruction, the model of the isolated dynamics of a single unit, a pair, a group of three units, etc. can be derived using proper identification techniques [34].
Our reconstruction technique works as follows.Suppose we have access to M measurements of the variables x 1 (t), . . ., x N (t) at times t respectively equal to h∆T with h = 1, ..., M , and ∆T being a (constant) sampling interval.The discretized version of Eqs.(1) reads: where i = 1, . . ., N , h = 1, . . ., M − 1, and x i (h) is a short notation for x i (h∆T ).Now, let y i (h
i,1 (M ) . . .σ 1 g (1) i,i+1 (M ) . . .σ 1 g i,N (M ) σ 2 g (2) where we introduced the following short notation: g i,j1,...,j d (h) := g (d) (x i (h), x j1 (h), . . ., x j d (h)).For each node i we need to identify terms, corresponding to the entries of A i .Therefore, Φ i ∈ R M ×H .Solving Eq. (3) for the unknowns A i allows one to reconstruct all interactions of node i, such that the whole structural connectivity can be inferred by repeating the calculations for all nodes, i = 1, . . ., N .Notice that Eq. ( 3) maps the problem of the reconstruction of the higher-order interactions into that of solving a system of algebraic equations in the unknown variables given by the H entries of A i .
Notice that, when M < H, the system of equations ( 3) is underdetermined and multiple solutions may exist [14].Conversely when M > H, the system of equations ( 3) is overdetermined and can be solved using the method of least squares [14].
We will now show that our approach is able to successfully reconstruct the full set of interactions at any order in the case of two completely different dynamics, namely microbial ecosystems and systems of coupled chaotic oscillators.The considered case studies will also demonstrate that the method works for the reconstruction of both hypergraphs and simplicial complexes, and no matter whether the underlying structure is undirected or directed, unweighted or weighted.

III. LOTKA-VOLTERRA DYNAMICS
In our first application we focus on the dynamics of microbial ecosystems.These consist of species that may engage in diverse relationships, either cooperative, such as the transfer of complementary metabolites, or antagonistic, such as competition for a resource [36].The validation of community-wide interactions in microbial communities is a far from trivial problem, faced both with experimental approaches [37] and through the use of mathematical modeling [38].The problem is further complicated by potential higher-order interactions, which play a role in stabilizing diverse ecological communities and maintaining species coexistence [1,39,40].Here, we model a microbial ecosystem of N species as a hypergraph of N coupled Lotka-Volterra equations [41] including both pairwise and three-body interactions: with i = 1, . . ., N .The variable x i represents the abundance of species i.The local dynamics of x i is governed by the logistic function f i (x i ) = r i x i 1 − 1 ki x i where r i and k i are the growth rate and the carrying capacity.The pairwise interactions between species are encoded in the coefficients of the N × N weighted matrix A (1) = {a ijk }.Note that Eqs. ( 5) are in the form of Eqs.(1) with g (1) (x i , x j ) = x i x j and g (2) (x i , x j , x k ) = x i x j x k .As an example, we consider the system of N = 7 species with four cooperative (a (1) ij > 0) and four antagonistic (a (1) ij < 0) pairwise interactions, studied in Ref. [38] and shown in Fig. 1(a) with blue and red arrows respectively.In addition to these pairwise interactions, we have included two antagonistic three-species interactions, shown as double arrows in the hypergraph in Fig. 1(a).These respectively correspond to a contribution to the dynamics of x 2 given by a (2) 237 x 2 x 3 x 7 and one to x 4 given by a (2) 416 = −0.0016[1].The other system parameters have been set as in [38].Namely, growth rates r i for all species have been selected from a uniform distribution in the interval (0, 1), similarly the carrying capacities k i from a uniform distribution in the interval (1, 100), and the initial conditions in the interval (10,100).
Under these conditions, as shown by the time evolution of the variables x i (t), with i = 1, . . ., 7, reported in Fig. 1(b), the microbial ecosystem typically converges to a stable equilibrium point corresponding to the coexistence of six species over seven.To feed our reconstruction algorithm, we sampled the time evolution of {x i (t)} i=1,...,7 at M regular intervals of size ∆T = 0.01, and we then used these values to calculate Y i and Φ i from Eq. (4).At this point, optimization via the method of least squares produces, for each i = 1, . . ., 7, a vector Âi that is solution of Eq. (3).To quantify the accuracy of the reconstruction of the interactions at any order, we compare the estimation Âi with the true values of the couplings, A i , for each i, evaluating the reconstruction error Fig. 1(c) shows E as a function of M/H.Different values of M/H have been obtained by changing the number of measurements M , while the number of unknown interactions to reconstruct H, given the symmetry of the interaction terms in Eq. ( 5), is fixed to The results indicate that our approach correctly reconstructs both pairwise and three-body interactions of the hypergraph, as the error vanishes when M/H > 1, i.e. when the system in Eq. (3) becomes overdetermined.This is evident also from Fig. 1(d) reporting the true values of the weights associated to non-zero interactions in the microbial ecosystem along with the values reconstructed for the case M/H = 0.95, when the system in Eq. ( 3) is underdetermined, and for the case M/H = 1.23, when the system is overdetermined.

IV. COUPLED R ÖSSLER OSCILLATORS
As a second case study we analyze the following system of Rössler oscillators coupled with pairwise and three-body interactions: where g (1) (x i , x j ) = x j − x i and g (2) (x i , x j , x k ) = x 2 j x k − x 3 i .As for the underlying topology of the interactions we consider simplicial complexes constructed as follows.We start from the so-called Zachary karate club, which is a system originally described in terms of a graph with N = 34 nodes and 78 links [42].Since the links form 45 triangles, we can represent the system as a simplicial complex by turning a fraction δ of the triangles into two-dimensional simplices [4].By considering different values of δ, we can then tune the percentage of the nodes forming a triangle which are effectively involved in a three-body interaction rather than in three, separate, pairwise interactions only.An example of the symplicial complex obtained for δ = 0.5 is shown in Fig. 2(a).Given that the interactions to reconstruct are, in this case, unweighted, we can now use three different methods to solve Eq. ( 3).In the first method, we consider ordinary least square minimization, namely we minimize the least square error between Y i and Φ i A i , as for the microbial ecosystem example.The green curves reported in Fig. 2(b) for δ = 0.5 and in Fig. 2(c) for δ = 1 show that the method correctly reconstructs the simplicial complex when M/H ≥ 1.In the second method, we consider minimization of the least square error under the additional constraint that the elements of A i are nonnegative: min Ai≥0 Y i − Φ i A i 2 .The orange curves in Fig. 2(b) and (c) indicate that, including such a priori information on the nature of the interactions in the optimization problem, expands the range of values of M/H for which reconstruction is possible.Lastly, we extend the signal lasso method [23] to deal with higher-order interactions.Namely, we consider the following optimization: min( Y i − Φ i A i 2 + α A i 1 + β A i − 1 H 1 ), where the penalty function includes, together with the 2-norm of the difference between Y i and Φ i A i , a regularization term to enforce sparsity of the solution, and another term to shrink the estimates to one.As indicated by the blue curves in Fig. 2(b) and (c), this method provides successful reconstruction in a larger range of values of M/H than the ordinary method of least squares.In conclusion, by using the last two methods, we have found that a smaller sample size is sufficient to fully reconstruct the structure of the simplicial complex.
V. CONCLUSIONS Summing up, in this work we have introduced an optimization-based framework to fully reconstruct the high-order structural connectivity of a complex system.Our approach can be useful to understand and predict the behavior of microbial ecosystems and coupled nonlinear oscillators, and we hope that it can shed new light on a variety of physical phenomena where higher-order interactions have a fundamental role.

FIG. 1 .
FIG. 1. Reconstructing higher-order interactions in a microbial ecosystem.(a) The underlying weighted hypergraph of a Lotka-Volterra system with N = 7 species and 2-and 3-body interactions, which we want to reconstruct from (b) the time evolution of the seven species abundance xi(t).(c) Quality of the reconstruction is measured by reporting the error E as a function of the ratio between the length M of the trajectories and the number H of interactions to reconstruct.(d) The reconstructed weights (color coded) for each 2-and 3-body interaction are compared to the true weights, for the two cases M/H = 0.95 ( ) and M/H = 1.23 ( ), also indicated in panel (c).

FIG. 2 .
FIG. 2. Testing the reconstruction method on a system of N = 34 coupled Rössler oscillators.(a) The underlying simplicial complex is obtained turning a tunable fraction δ of the triangles of the Zachary karate club network into 2-simplices (case δ = 0.5 is shown).(b-c) Reconstruction error E as a function of M/H for δ = 0.5 and δ = 1.0.

( 1 )
ij }, while the three-body interactions in the coefficients of the N × N × N weighted tensor A(2) = {a(2)