Reducibility of higher-order networks from dynamics

Domain: physics
Published: 2026-01-15
Editorial outcome: major_revision

Abstract
Empirical complex systems can be characterized not only by pairwise interactions, but also by higher-order (group) interactions influencing collective phenomena, from metabolic reactions to epidemics. Nevertheless, higher-order networks’ apparent superior descriptive power—compared to classical pairwise networks—comes with a much increased model complexity and computational cost, challenging their application. Consequently, it is of paramount importance to establish a quantitative method to determine when such a modeling framework is advantageous with respect to pairwise models, and to which extent it provides a valuable description of empirical systems. Here, we propose an information-theoretic framework, accounting for how structures affect diffusion behaviors, quantifying the entropic cost and distinguishability of higher-order interactions to assess their irreducibility to lower-order structures while preserving relevant functional information. Empirical analyses indicate that some systems retain essential higher-order structure, whereas in some technological and biological networks it collapses to pairwise interactions. With controlled randomization procedures, we investigate the role of nestedness and degree heterogeneity in this reducibility process. Our findings contribute to ongoing efforts to minimize the dimensionality of models for complex systems.

Introduction
Many complex systems exhibit an interconnected structure that can be encoded by pairwise interactions between their constituents. Such pairwise interactions have been used to model biological, social and technological systems, providing a powerful descriptive and predictive framework1–6. Recently, the analysis of higher-order structural patterns and dynamical behaviors attracted the attention of the research community as a powerful framework to model group interactions7–10, with applications ranging from neuroscience11,12 to ecology13 and social sciences14,15, highlighting the emergence of novel phenomena and non-trivial collective behavior16–21.

Higher-order networks encode more information than pairwise interactions: for example, metabolic reactions are more realistically described by group interactions between any number of reagents and reactants, capturing information that would be lost by considering the union of pairwise interactions instead. However, this modeling flexibility comes at a cost: new data needs to be adequately recorded and stored as group interactions instead of pairwise ones, and new analytical22–25 and computational26–28 tools need to be developed. Moreover, the complexity and computational cost of these tools increase exponentially as larger group interactions are considered. It is therefore crucial to understand under which conditions higher-order representations need to be favored over classical pairwise ones, and whether it is possible to devise a quantitative procedure to determine which representation provides a suitable description of an empirical system.

A research line approaches higher-order interactions from the perspective of “behavior”, seeking to characterize statistical dependencies in multivariate time series29–33. Within this framework, information-theoretic methods have shown that such dependencies can be decomposed into contributions from different interaction orders34. At the same time, initial studies have revealed that the relationship between higher-order behaviors and mechanisms, i.e., structure and dynamics, is subtle and often non-trivial10,35. How to reduce a higher-order structure (and the dynamics it supports) given the dynamics on top of it, remains an open problem.

A similar challenge was faced nearly a decade ago, when the advent of temporal and categorical data36–38 allowed multilayer representations of complex networks39. For these representations, an information-theoretic approach was used to show that not all layers, or types of interaction, are equally informative: Some information can be discarded or aggregated to reduce the overall complexity of the model without sacrificing its descriptive power40,41. Although multilayer networks are different from higher-order networks, this approach, formally based on the density matrix formalism42, provides a good candidate for the present case. Indeed, the idea is similar in spirit to the widely used information compression algorithms adopted in computer science: by exploiting the regularities in the data, one can build a compressed representation that optimizes the number of bits—i.e., the description length—needed to describe the data with a model and thus to encode the model itself. Broadly speaking, this minimum description length principle, grounded on Bayesian inference, allows one to maximize model accuracy while minimizing model complexity by means of a principled regularization term43–45. Similar approaches have also been used to coarse-grain complex and multiplex networks46–48.

Here, we build on this long-standing research line and propose an approach to compress systems with higher-order interactions, while accounting for the ability to describe complex data and the cost due to the complexity of the model. Specifically, we determine a suitable order up to which interactions need to be considered to obtain a reduced and functionally plausible representation of the system. Larger orders above this threshold can be safely discarded. Formally, we do so by generalizing the concept of network density matrix42,49 to account for higher-order diffusive dynamics with the multiorder Laplacian17,50 and deriving a meaningful rescaling of the diffusion time as a function of the order of interaction. Then, we calculate an information-theoretical cost function which depends on the largest order considered. By minimizing this cost function, we find the optimal compression of the data which, in turn, corresponds to the most parsimonious functional description of the system, at a given diffusion time. In the following, we refer to this procedure as functional reduction. A higher-order network is fully reducible—to a pairwise network—if the optimal order is 1, while its reducibility decreases for increasing optimal orders. To this aim, we consider two alternative, yet complementary, perspectives: the first, based on the Kullback-Leibler divergence DKL(ρ||ρ0), quantifies the entropic cost of building a system represented by a density matrix ρ from a collection of isolated nodes represented by a density matrix ρ0; the second, based on the Jensen-Shannon divergence DJS(ρ||ρ0), measures how distinguishable the system is from such a collection. These two measures are complementary alternatives, each with distinct advantages and limitations in capturing reducibility, as they emphasize different informational aspects of the system. The choice between them depends on the specific application and context. For simplicity, in the following we present results based on the former, and refer to the Supplementary Information for the analysis using the latter.

We analytically derive the cost function for a rotationally invariant structure, revealing how reducibility depends on Laplacian eigenvalues and diffusion time. Then, we demonstrate the validity of our method by performing an extensive analysis of synthetic networks and investigate the functional reducibility of a broad spectrum of real-world higher-order systems. Our analysis reveals (i) that datasets from different domains are differently reducible, (ii) that no simple structural metric alone can explain reducibility, and (iii) the role of degree heterogeneity in reducibility by using randomization procedures.

The advantage of this framework is that it provides a bridge between network analysis and information theory by means of a formalism that is largely inspired by quantum statistical physics, which has found a variety of applications from systems biology51 to neuroscience52, shedding light on fundamental mechanisms such as the emergence of network sparsity53 and the renormalization group54.

Results
Density matrix for higher-order networks
The flow of information between nodes in a (pairwise) network can be modeled by different dynamical processes. However, in the absence of specific knowledge or hypotheses about the underlying dynamics, diffusion represents the maximum-caliber estimate and thus provides a principled modeling choice53. Diffusion on networks can be described by means of the propagator e−τL, where L the combinatorial Laplacian and τ is the diffusion time. In particular, the information flow from node i to node j is described by the component (e−τL)ij. Network states can then be encoded by a density matrix42,49 defined as

ρτ = e−τL / Z, (1)

where the partition function Z=Tr(e−τL) ensures a unit trace for this operator. The influence through the network guided by the propagator is often described as information flow. However, we note that this should not be conflated with information flow in the information-theoretic sense, as quantified by entropy, mutual information, or transfer entropy. Here, the term is rather used in the broader network science convention to indicate the spread of dynamical effects through the system.

Here, we generalize density matrices to higher-order networks. The most general formalism to encode higher-order networks is that of hypergraphs8. A hypergraph is defined by a set of nodes and a set of hyperedges that represent the interactions between any number of those nodes. A hyperedge is said to be of order d if it involves d+1 nodes: a 0-hyperedge is a node, a 1-hyperedge is a 2-node interaction, a 2-hyperedge is a 3-node interaction, and so on. Simplicial complexes are a special case of hypergraph that is also commonly used: they additionally require that each subset of each hyperedge is included in the hypergraph, too. In a hypergraph H with maximum order dmax, diffusion between nodes through hyperedges of order up to D≤dmax can be described by the rank–2 projected operator known as multiorder Laplacian50,55

L[D] = L(D,mul) = Σd=1D γd L(d) / 〈K(d)〉, (2)

which is a weighted sum of the Laplacians L(d) at each order up to order D. At each order, the weight is defined by a real coefficient γd (which we set to 1 for simplicity) and the averaged generalized degrees 〈K(d)〉. Each d-order Laplacian is defined by L(d)ij = K(d)i δij − A(d)ij, in terms of the generalized degrees K(d) and adjacency matrix A(d) of order d. The matrix L[D] satisfies all the properties expected from a Laplacian: it is positive semidefinite and its rows (columns) sum to zero (see Methods Section I for details). In general, a hypergraph does not need to have hyperedges at every order below D, unless it is a simplicial complex. If there is no hyperedge of order d, both the Laplacian and the averaged degree in Eq. (2) vanish, and the result is undefined. In those cases, the sum thus needs to be taken over all orders below D that exist: D = {d≤D: 〈K(d)〉>0}.

The multiorder density matrix of hypergraph H, up to order d, is defined as

ρτ[d] = e−τL[d] / Z, (3)

with the partition function Z=Tr(e−τL[d]). Just like its pairwise analog in Eq. (1), this operator satisfies all the expected properties of a density matrix: it is positive definite, and its eigenvalues sum up to one. Importantly, the diffusion time τ allows playing the role of a topological scale: small values of τ allow information to diffuse only to neighboring nodes, probing only short-scale structures. Larger values of τ, instead, allow the diffusion to reach more remote parts of the hypergraph and describe large-scale structures. In this context, the meaning of “small” and “large” depends on the network structure, and can be estimated with respect to the magnitude of the largest (1/λmax) and smallest (non-zero) (1/λmin) eigenvalues of the Laplacian matrices, respectively. To meaningfully compare the density matrix in Eq. (3) at different orders d, we will need to rescale the chosen τ at each order—the derivation and necessity of the rescaling are detailed below in Section C.

Quantifying the reducibility of a hypergraph
We approach the reducibility of a hypergraph as a problem of model selection. We formulate that problem as follows: given a hypergraph H with maximum order dmax, is H an optimal representation of itself, or is considering only its hyperedges up to a given order d<dmax sufficient? Formally, we treat the density matrix ρ[dmax] as data and ρ[d] as a model of the data. Formulated in this way, we need to find the “optimal” model of the data, that is, to evaluate the optimal order dopt. To this aim, it is desirable to frame the optimization problem in Bayesian terms, following the minimum description length principle, akin to Occam’s razor. That framework quantifies the suitability of a model through its likelihood and its complexity through the prior. However, one important limitation should be noted in our case: defining a prior in terms of the density matrix is non-trivial, preventing the definition of a unique and self-consistent complexity term as in the original minimum description length principle. To address this, we introduce suitable—though not necessarily optimal—cost functions. Another important aspect is that the likelihood function provides only a proxy description of the observed network, expressed through the diffusive process unfolding on it. This requires interpreting our results in terms of the interplay between structure and dynamics, rather than structure alone.

The steps of the methods are illustrated in Fig. 1: (a) start from an original hypergraph, (b) use higher-order diffusion at a selected diffusion time τ, then (c) calculate the optimal largest order dopt by minimizing the cost function, and (c) reduce the original hypergraph to a suitable hypergraph, that is, one with orders only up to dopt without losing functional information.

As illustrated in Fig. 2, we define the cost function L as the sum of the information loss—the opposite of the model accuracy—of the model, measured by a suitably generalized Kullback–Leibler divergence DKL, and the model complexity C (see “Methods” Section K for details):

L(ρτ[dmax]|ρτ[d0]) = DKL(ρτ[dmax]|ρτ[d0]) + C(ρτ[d0]), (4)

where τ0 is a necessary order-dependent rescaling of τ, as will be explained in Section C.

Note that by definition, there is no information loss when considering all possible orders, i.e., DKL(ρτ[dmax]|ρτ[dmax]) = 0. It is worth remarking that this cost function deviates from the standard Bayesian description length for probability distribution, because it is not trivial to define the corresponding prior term in terms of operators (see “Methods”). To partially overcome this limitation, we use again a Kullback-Leibler divergence to approximate the model complexity as C(ρτ[d]) = DKL(ρτ[d]|ρiso) = logN − Sτ[d], where Sτ[d] = −Tr(ρτ[d]logρτ[d]) is the von Neumann entropy of the hypergraph up to order d, and ρiso = I/N is the density matrix of N isolated nodes (see “Methods” Section L). The complexity can take values in C(ρτ[d]) ∈ [0, Siso], reaching its minimum for hypergraphs with maximal entropy (e.g., isolated nodes or hypergraphs with high regularity), and its maximum for hypergraphs with minimal entropy (e.g., hypergraphs with more complex structures).

Another definition, grounded on quantifying distinguishability, is discussed in Section M and described in S.I. Sec. V.

In the following, we will refer to optimal order dopt, at diffusion time τ, as the most suitable order that minimizes our cost function:

dopt(τ) = argmin d L(ρτ[dmax]|ρτ[d0]) : (5)

It is also worth remarking that the optimal order should be interpreted in terms of the diffusion dynamics on the network, not directly its structure. With this prescription, the likelihood term—encoded by DKL—in the cost function retains all the features of a likelihood function42 and can be suitably used for our purpose.

Finally, we define the reducibility of the hypergraph at diffusion time τ as

χτ(H) = (dmax − dopt) / (dmax − 1), (6)

which measures the ratio between the number of orders to reduce and the maximum number of orders to reduce, dmax − 1. By construction, χτ(H)=0 for a hypergraph that is not reducible at all, i.e., dopt(τ) = dmax, while χτ(H)=1 for a hypergraph that is maximally reducible, that is, it can be optimally reduced to its pairwise interactions, dopt(τ) = 1. The dependence on τ makes this a multiscale method: each diffusion time serves as a different lens to probe the system. There is typically a range of informative τ values determined by the eigenvalues of the multiorder Laplacian, as we will see in Section D.

Rescaling the diffusion time τ at each order
As mentioned above, a given diffusion time τ may be large for some networks but low for others, which is a challenge when the aim is to compare networks and orders. To overcome this problem, we rescale τ to ensure an appropriate diffusion time for each structure, and we exploit examples of hypergraphs with certain regularities to define a baseline and characterize the scaling relation.

Specifically, there is a class of hypergraphs for which Laplacians of different orders are proportional, i.e., L(d) ∝ L(d0). This occurs in very regular structures, including complete hypergraphs and some simplicial complex lattices (see “Methods” Section N). In this case, since the Laplacian matrices govern the flow of information, all orders and all their combinations encode the same functional information. Consequently, we expect these structures to be functionally invariant under reduction—i.e., the cost function should not change as one reduces the hypergraph.

However, without rescaling, the summation of Laplacian matrices according to Eq. (2) simply strengthens the flow pathways in the original hypergraph, making it different from its reduced versions. To correct for this effect, we rescale τ to allow meaningful comparisons between a hypergraph at different orders.

Since the density matrix only depends on the product of τL(d), this can be achieved by selecting a value of diffusion time τ, and then rescaling it to obtain a new τ0(d) at each order like so:

τ0(d) = (dmax / d) τ, (7)

which we denote τ0 ≡ τ0(d) to avoid cluttering the notation. This rescaling ensures that the multiorder density matrices are all equal, ρτ0[d] = ρτ[dmax] ∀d, in the special case of hypergraphs with proportional Laplacian matrices, and gives a flat cost function L(ρτ[dmax]|ρτ[d0]) = C(ρτ[dmax]) for all d in those extreme structures (S.I. Sec. I).

For illustration purposes, we set τ = τshort = 1/λmax in numerical experiments, unless otherwise stated, where λmax is the largest eigenvalue of the multiorder Laplacian of the original hypergraph L[dmax]. This corresponds to choosing a small diffusion time to explore more local structures of the hypergraph.

Physically, this rescaling simply means that we adapt the diffusion time at each order so that the diffusion process on the hypergraph always reaches the same accumulated strengths. We use this rescaling of τ at each order in all hypergraphs. With this prescription, the likelihood term—encoded by DKL—in the cost function retains all the features of a likelihood function42 and can be suitably used for our purpose.

Analytically tractable case: the hyperring
For arbitrary hypergraphs, the cross-entropy part of the information loss is hard to derive analytically. Here, we consider a case that is analytically tractable because of its rotational symmetry: a hyperring56 with interactions up to order 2.

We define the hyperring as N nodes on a ring with pairwise interactions with first neighbors, (i,i+1)modN, ∀i and triplet interactions with first neighbors too (i−1,i,i+1)modN, ∀i. By construction, this hyperring is a simplicial complex. Because of the rotational symmetry, the Laplacian and density matrices are circulant matrices, which commute, and for which the eigenvalue spectrum is known and expressed in the Fourier basis56,57. Specifically, the multiorder Laplacian spectra are (see full derivation in Supp. Info. Sec. II), up to order 2:

λ[1]k = 1 − cos(2πk/N), (9)

λ[2]k = 2(1 − cos(2πk/N) − cos(4πk/N)), (10)

with λ[d]k = 0 and λ[d]k = λ[d]N−k for 0 < k < ⌊(N−1)/2⌋ and d=1,2 (Fig. 3a). We also know that the spectra of the corresponding density matrices can be expressed in terms of these as μ[d]k = e−τ0λ[d]k / Z[d], which depend on τ (Fig. 3b). These eigenvalues sum to 1 and can be seen as probabilities, or weights associated with each diffusion mode. For short diffusion times, these weights are distributed uniformly (blue curves), but at longer diffusion times, the fast modes decay, and only the slow modes remain (green curves).

We can now explicitly write the information loss: indeed, it reduces to a KL divergence over the distributions of eigenvalues μ[d]k because the matrices are circulant and hence commute. Eventually, we can write

DKL(ρτ[dmax]|ρτ[d0]) = τ(dmax−1) Eρ[dmax](λ[d]k) − Eρ[dmax](λb[d]k) + log(Z[d]/Z[dmax]), (11)

where we denote

Eρ[dmax](λ[d]k) = ΣN−1k=0 μ[dmax]k λ[d]k, (12)

the expectation of a given spectrum weighted by the diffusion modes under the full system diffusion, and

λb[d]k = δd,dmax + (1 − δd,dmax) Eλ[d](δ), (13)

the eigenvalues of the hypergraph composed by all truncated orders. Eq. (11) is interpretable and has an overall structure similar to those derived for pairwise networks42. The first two terms are proportional to τ and account for “spectral distortion”: (1) contributed by orders up to d, and (2) not contributed by orders d<δ≤dmax not present in the truncated system. The (dmax−1) pre-factor accounts for the rescaling of τ and vanishes when d=dmax. The third term is a standard log-correction of the spread of the two density matrix spectra. This information loss vanishes in the limits τ → 0 and τ → ∞ (in the former, the expected values vanish because diffusion peaks in the uniform mode k = 0, the only one not decaying). This indicates that these extreme diffusion times are not informative scales at which to probe the system: at τ = 0 no diffusion has occurred, whereas at τ = ∞ all diffusions have reached the uniform equilibrium. At intermediate τ, the information loss vanishes for d=dmax, as expected.

We can now derive the complexity in terms of those same eigenvalues:

C(ρτ[d0]) = log Z[d] − τ(dmax−1) Eρ[dmax](λ[d]k). (14)

Note that this expression of the complexity is valid for any hypergraph, as the symmetry of the hyperring was not needed to derive it. The complexity tends to zero, C(ρτ[d0]) → 0 for τ → 0, and C(ρτ[d0]) → logN for τ → ∞.

The cost function is simply the sum of Eqs. (11) and (14). In order to obtain the optimal order as a function of the diffusion time τ, we can simply compute the difference between the cost function at orders 2 and 1 (Fig. 3c): the hyperring is irreducible at short diffusion times (dopt=1) but reducible at intermediate and long diffusion times (dopt=2). The sign of the result indicates whether the optimal order is 1 or 2. Note again that extremely short τ → 0 or long τ → ∞ diffusion times are not informative scales at which to probe the system.

These analytical expressions are interpretable and reveal how reducibility is directly linked to the structure and diffusion dynamics on top of it via the diffusion time, and the eigenvalues of the multiorder Laplacian. Explicitly, the terms Eρ[dmax](λ[d]k), Eρ[dmax](λb[d]k) quantify the spectral distortion between the truncated and the full hypergraph. At τ ~ τshort, all modes contribute nearly equally, so reducibility depends on whether higher-order edges introduce genuinely new local features into the spectrum. At τ ~ τlong, only the slow modes dominate, making reducibility hinge on whether higher orders reshape the global connectivity patterns. In both regimes, the trade-off with model complexity determines whether the gain in descriptive power outweighs the penalty of keeping additional orders. Hence, whenever adding higher-order edges does not substantially alter the relevant spectral content at the probed scale, the system will be deemed reducible. This interpretation also explains why the hyperring appears reducible at most (not-too-small) diffusion times: the added higher-order edges do not substantially alter global patterns of diffusion, because of the rotational symmetry and nestedness—strong structural constraints. From this point of view, the hyperring provides a conceptual bridge between the complete hypergraph, and the more complex, and less symmetric, structures that we investigate numerically in the next sections.

Random structures
We now investigate the reducibility of two types of heterogeneous random structures: random hypergraphs and random simplicial complexes. The models are described above. In both cases, we set N=100, pd=50/Nd and dmax=4.

Figure 4 shows the cost function considering orders from 1 to 4, for 100 realizations of each type of random structure, at τ = τshort. In the case of random simplicial complexes (Fig. 4a), dopt=3, reflecting some reducibility χ = 1/3. Instead, in random hypergraphs (Fig. 4b), the optimal order is the maximum, dopt=4=dmax. This mean that those random hypergraphs, at this scale, are not reducible, i.e., χ=0.

The crucial difference between these two cases is that the hyperedges between different orders are correlated (nested) in random simplicial complexes but uncorrelated in random hypergraphs. Indeed, simplicial complexes need to satisfy downward inclusion: if a d-hyperedge is present, all (lower order) hyperedges that are subsets of its nodes need to be included too. More information about the differences between the two and their impact on dynamics can be found in ref. 17. Nestedness can be measured by the simplicial fraction defined in ref. 58: it is 1 (maximal) in simplicial complexes but closer to 0 (minimal) in random hypergraphs.

To test the effect of nestedness on reducibility, we start from a random simplicial complex (Fig. 4c) and gradually change it into a random hypergraph by shuffling its hyperedges. Specifically, we change each hyperedge into another non-existing hyperedge with probability pshuffle, preserving the number of hyperedges. For pshuffle = 1, the result is a random hypergraph (Fig. 4e). Figure 4c–e show how decreasing nestedness also decreases reducibility.

In simplicial complexes, higher-order hyperedges are strongly nested within lower-order ones. Therefore, higher-order edges that overlap with lower-order ones add little new spectral content, yielding reducibility. By contrast, in random hypergraphs or when nestedness is destroyed through shuffling, additional hyperedges substantially reshape the spectrum by introducing genuinely new local connections, thereby lowering reducibility.

As mentioned above, τ acts as a topological scale. Diffusion processes with different values of τ are thus expected to “see” different structures, which may result in different optimal order and reduced hypergraph. To illustrate this, we compute the cost function curves for the random simplicial complex case, for five values of τ evenly spaced on a logarithmic scale, and where the second and fourth are 1/λmax and 1/λmin, respectively. Supplementary Fig. S2 shows that as the diffusion time increases (larger τ), the cost function (i) increases overall and (ii) the reducibility decreases to χ=0. Supplementary Fig. S4 shows the same results for three values of the hyperedge density, where we can see that for higher densities, the cost function curve becomes much flatter with no clear minimum, suggesting a functional similarity between the large-scale structure at all orders.

We similarly tested the effect of the sole density on the reducibility (Supplementary Figs. S3 and S4). In general, the overall density coefficient does not appear to significantly affect the reducibility.

At long diffusion times, reducibility is generally low because the expectations weight only the slow modes, which capture global connectivity patterns. In this regime, higher-order edges have a stronger effect: unless they are tightly nested within lower-order ones, they rewire long-range connections and shift the smallest eigenvalues, making the truncated spectra diverge from the full one. As a result, simplicial complexes retain some reducibility thanks to their nested structure, but random hypergraphs or shuffled ensembles remain irreducible at long times.

Real-world hypergraphs display diverse reducibility
We now investigate how reducible real-world hypergraphs are by considering 60 empirical hypergraph datasets from 9 categories: coauthorship, face-to-face contacts, drugs, ecology, email, metabolism, politics, online forum tags, and theater.

First, we observe a great variety in the reducibility values (Fig. 5) and the associated cost function curves (Supplementary Figs. S5 and S6) for τ = τshort. Figure 5a–c shows three examples with curves of very different shapes and with different reducibility. All reducibility values are shown in Fig. 5d by category, and their overall distribution cover the full possible range with many values close to 0 (Fig. 5e). Second, we note that datasets from different categories seem to be distinctively reducible: none of the coauthorship, email, theater, and ecology datasets are very reducible (χ≈0), whereas the contact, drugs, and metabolism datasets show a range of reducibility values, and the tags datasets appear fully reducible (χ = 1). The variability between the drugs datasets can be understood because the nodes and hyperedges represent different concepts in each case (see “Methods” Section O). A Kruskal-Wallis H-test indicated that the category of the dataset has a statistically significant effect on the reducibility, p<10−4—after excluding the categories with too few (≤5) datasets. A description of all datasets can be found in “Methods” Section O, and their summary statistics and reducibility can be found in Table 1.

For completeness, we also show the reducibility values for a long τ = τlong in Supplementary Fig. S7. We observe that in most cases, the reducibility drops to zero in that case, as in the synthetic cases. This can be understood as follows: At long diffusion times, the dynamics is governed by the slow Laplacian modes (small eigenvalues), which capture global connectivity patterns. Adding higher-order edges almost always alters these global patterns—by creating shortcuts, bridging distant parts, or biasing flows across large subsets of nodes—so the truncated spectra diverge strongly from the full one. As a result, information loss dominates and the system becomes irreducible.

Reducibility is not explained by any one simple metric
To better understand the relationship between reducibility and the structure of the hypergraphs, we computed the Pearson correlation coefficient between the reducibility and several typical structural metrics.

In Fig. 6, we show the results for six metrics: density (M/N, where M is the number of hyperedges and N the number of nodes), the logarithm of the maximum generalized degree, the logarithm of the spectral radius (that is, the maximum of the absolute eigenvalues of the generalized adjacency matrix), nestedness (as measured by the simplicial fraction58), the cross-order degree correlation (between orders 1 and 2), and the heterogeneity of the generalized degree (measured as (Kmax − 〈K〉)/Kmax). Note that although all metrics use quantities generalized for higher-order networks, such as degree and density, nestedness and cross-order degree correlation are strictly higher-order metrics: they measure how the structure at different orders relate. For each metric, we report the correlation coefficient r and the associated p-value.

Results indicate a significant (p < 0.05) positive correlation between reducibility and each of these metrics except for the cross-order degree correlation. However, note that in all cases, for a given value of a structural metric, reducibility still displays a range of values. For example, reducibility is correlated to density (r=0.51, p=3.5×10−5) but density alone cannot explain reducibility: datasets with low density exhibit reducibility values between 0 and 1. Note that this might in part be because most of these metrics only consider the structure aggregated across orders, whereas reducibility considers how the structure and dynamics change as more orders are considered. One exception of this is the nestedness, which accounts for the relation between orders (and cross-order degree correlation, but here only computed between orders 1 and 2 for practical reasons). Hence, if a hypergraph has high nestedness, it gives us information about how the structure changes as larger orders are added.

For completeness, we report the correlation between reducibility and several additional metrics in Supplementary Fig. S8a–f and the correlations between those metrics in Supplementary Fig. S8g.

Linking reducibility to degree heterogeneity by randomizing
We start by focusing on nestedness as a first metric to help us better understand reducibility in empirical datasets. Nestedness describes how often hyperedges of lower orders are subsets of hyperedges of higher orders. Simplicial complexes are maximally nested because they satisfy closure downwards, whereas random hypergraphs are minimally nested because there is no correlation between its hyperedges at different orders. To quantify nestedness, we use the simplicial fraction as above.

There are at least three reasons to start with nestedness: (i) nestedness is a form of structural redundancy, (ii) in synthetic hypergraphs, we observed that decreasing nestedness in a controlled manner decreased reducibility (Fig. 4), and (iii) in empirical datasets, we observed a positive correlation between nestedness and reducibility. As noted in ref. 58, many empirical datasets are not very nested. This is also the case in the datasets considered here, with the only nested ones belonging to the “contact” category. For this reason, we now focus on that category, which also displays a range of reducibility values (Fig. 5).

To investigate the link between reducibility and nestedness, we use three randomization strategies (from least to most structure-destructive): the configuration model59, a node swap procedure17, and a simple shuffling of hyperedges17 (see “Methods” Section I for details).

In the configuration model, we randomize each order of interaction independently. This randomizes hyperedges while preserving the degrees of the nodes at each order. In the node swap procedure, at each order d>1, we pair nodes (i,j) and swap their generalized degrees K(d)i ↔ K(d)j by swapping their hyperedge memberships at that order. This preserves the degree distributions at each order, but destroys the degree sequences (i.e., the assignment of degrees to specific nodes). In the shuffling, we simply shuffle all hyperedges. This randomizes hyperedges too, while only preserving the number of hyperedges—the degree sequences and distributions are destroyed. This shuffling was used in the synthetic random structure (Fig. 4c, d). Compared to the first two randomizing procedures, the shuffling is extreme and yields structures akin to random hypergraphs. Theoretically, we expect these three randomization procedures to destroy nestedness in the datasets since they destroy correlations between orders, and we are interested in seeing their effect on reducibility.

We show the results for 20 realizations of each randomization procedure in Fig. 7. First, we observe that all three procedures destroyed nestedness, as expected (Fig. 7b): original values (black) were mostly 0.7 (and one around 0.4), but nestedness is close to 0 in all structures randomized by the configuration model (green), the node swap (blue), and the shuffling (red). As anticipated, the shuffling also caused the reducibility to vanish (Fig. 7a), consistent with what we observed in synthetic structures (Fig. 4c, d). Interestingly, however, the configuration model and the node swap preserved the reducibility values most of the time. In other words, in hypergraphs with originally large nestedness and reducibility, the reducibility is preserved if the degree sequences and distributions are, too, but is destroyed if not. This indicates that it is the information contained in the degree distributions at each order, rather than the nestedness, that is impacting the reducibility of these systems.

To go further, we examined the cross-order degree correlation and the degree heterogeneity. The first measures how correlated the degrees of nodes at different orders are (Fig. 7c). It is a summary statistic of the degree sequences (but cannot be computed from the degree distributions, as it contains no information about the degree-node assignments). As a proxy, we simply look at the correlation between orders 1 and 2: DC=corr({K(1)i}, {K(2)i}). Second, the degree heterogeneity measures how varied the generalized degrees are (Fig. 7d). It is a summary statistic of the degree distributions (and sequences, but the distributions are sufficient). In very nested hypergraphs, “richer nodes get richer” and cross-order degree correlation is high17. As expected, the original “contact” datasets have high DC (close to 1) which is destroyed (close to 0) by the shuffling and the node swap. DC is only preserved by the configuration model, that is, by preserving the degree sequences, as expected. This indicates that the DC cannot explain the reducibility in this case, since the reducibility can be preserved when the DC is destroyed (node swap). Then we turn to the generalized degree heterogeneity: it is only destroyed by the shuffling, and thus preserved by preserving degree sequences (configuration model) or distributions (node swap). We thus observe that it is sufficient that the randomization preserves degree heterogeneity to preserve reducibility.

We also checked the effect of randomizations on other metrics: the higher-order degree heterogeneity ratio17, the spectral radius, and maximum higher-order degree (Supplementary Fig. S9). The maximum degree, another summary statistic of the degree distributions, shows a similar but less marked trend than the degree heterogeneity: it decreases but is not fully destroyed when the shuffling destroys the reducibility. The degree heterogeneity ratio and spectral radius do not exhibit any clear trend in this case.

In conclusion, the results suggest that the degree distributions and their summary statistic of degree heterogeneity (and, to a lesser extent, the maximum degree) are key to the reducibility of these “contact” datasets that are highly nested, whereas nestedness and cross-order degree correlation are not. From the perspective of spectral distortion, this suggests that heterogeneity—which we know is directly linked to the broadness of the Laplacian spectrum—correlates in these datasets with not much new spectral content as orders are added, at short times, thereby enabling lower-order truncations to approximate the full spectrum. Applying other suitable randomizations to other categories of datasets (with lower nestedness) could help disentangle other key ingredients to reducibility. Interestingly, past research showed that cross-order degree correlation affected dynamics on hypergraphs: low DC stabilizes synchronization17, and higher DC facilitates bistability contagion models59.

Methods
Multiorder Laplacian for hypergraphs
The Laplacian matrix L of a graph is defined as Lij = Kiδij − Aij, where Ki is the degree of node i, δij is the Kronecker delta, and Aij is the adjacency matrix. For hypergraphs, we can define d Laplacian matrices, one for each order of interaction. The Laplacian of order d can be defined as17,50

L(d)ij = K(d)i δij − (1/d) A(d)ij, (15)

where K(d)i is the degree of order d of node i, i.e., the number of d-hyperedges connected to node i, while A(d)ij is the adjacency matrix of order d, whose elements A(d)ij counts the number of d-hyperedges connected to nodes i and j. We can define the multiorder Laplacian up to an order D as17,50

L[D] = L(D,mul) = Σd=1D γd L(d) / 〈K(d)〉, (16)

where γd is a tuning parameter of interactions of order d, while 〈K(d)〉 is the averaged degree of order d. For simplicity, in this study, we always set γd = 1.

Note that in general, a hypergraph does not need to have hyperedges at every order below D, unless it is a simplicial complex. If there is no hyperedge of order d, both the Laplacian of order d and the averaged degree in Eq. (2) vanish and the result is undefined. In those cases, the sum thus needs to be taken over all orders below D that exist: D = {d≤D: 〈K(d)〉>0}.

Multiorder density matrix for hypergraphs
We define the multiorder density matrix of a hypergraph H, up to order d, as

ρτ[d] = e−τL[d] / Z, (17)

with the partition function Z=Tr(e−τL[d]) and the diffusion time τ. As its pairwise analog in Eq. (1), this operator satisfies all the expected properties of a density matrix: it is positive definite, and its eigenvalues sum up to one.

We know from ref. 42 how the eigenvalues of the density matrix λi(ρτ[d]) are related to those of the Laplacian λi(L[d]):

λi(ρτ[d]) = e−τλi(L[d]) / Z, (18)

where the partition function can also be expressed as Z = ΣNi=1 e−τλi(L[d]). This relation is useful to link our intuition and knowledge about the Laplacian spectrum to that of the density matrix.

The diffusion time τ plays the role of a topological scale parameter: small values of τ allow information to diffuse only to neighboring nodes, probing only short-scale structures. Larger values of τ, instead, allow the diffusion to reach more remote parts of the hypergraph and describe large-scale structures. In this context, the meaning of small and large depends on the network structure, and can be estimated with respect to the magnitude of the largest (1/λmax) and smallest (1/λmin) eigenvalues of the Laplacians, respectively.

Information loss as Kullback–Leibler divergence between two hypergraphs
The state of the original hypergraph H is stored in the multiorder density matrix

ρτ[dmax] = e−τL[dmax] / Z[dmax], (19)

and the state of the reduced hypergraph where we consider orders up to d is given by

ρτ[d] = e−τL[d] / Z[d]. (20)

To perform model selection and determine the optimal order dopt to represent the hypergraph, we treat ρτ[dmax] as data and ρτ[d] as a model of it. The first key aspect of a good model is that it must describe the data as accurately as possible. We can quantify the modeling error, or information loss, with the Kullback–Leibler (KL) entropy divergence between the data and the model, defined as

DKL(ρτ[dmax]|ρτ[d]) = −S[dmax] + S(ρτ[dmax]|ρτ[d]) ≥ 0, (21)

where S[dmax] = −Tr(ρτ[dmax]logρτ[dmax]) is the Von Neumann entropy of the hypergraph and S(ρτ[dmax]|ρτ[d]) = −Tr(ρτ[dmax]logρτ[d]) is the cross-entropy between the hypergraph and its reduced form.

The von Neumann entropy of a hypergraph can also be written as

Sτ[d] = −Tr(ρτ[d]logρτ[d]) = − Σi λi(ρτ[d])logλi(ρτ[d]), (22)

where {λi(ρτ[d])} are the eigenvalues of the density matrix, and 0log0=0 by convention. The von Neumann entropy S[d] is zero if only one eigenvalue is non-zero (“pure state”, in the language of quantum mechanics), and maximal equal to logN if all eigenvalues are equal (“maximally mixed state”, in the language of quantum mechanics). This maximal entropy is realized in a trivial hypergraph with N isolated nodes42: all eigenvalues of the Laplacian are zero, so Z = N, and λi(ρτ[d]) = 1/N ∀i according to Eq. (18), yielding Sτ[d] = logN. Another extreme case is that of the complete hypergraph. In that case, one eigenvalue is λ1(Lτ[d]) = 0 and the other N − 1 are λi(L[d]) = dN/(N−1) (see Eq. (28) below). The expression for Sτ[d] can then be derived analytically and, consistent with the case of the complete pairwise network42, it has minimal entropy Sτ[d] = 0 in the limit τ → +∞ and maximal entropy Sτ[d] = logN in the limit τ → 0.

The information loss DKL takes positive values proportionally to the inaccuracy of the model, and reaches zero at d=dmax—as the data is the most accurate model of itself. It is to avoid this overfitting that a regularization term, accounting for model complexity, is needed.

Model complexity
It is compulsory that high model accuracy must be balanced with low model complexity.

Within the Bayesian framework underlying the minimum description length principle, model complexity arises from the prior term. However, as discussed above, defining a valid prior in terms of the density matrix is nontrivial. Accordingly, we adopt a more heuristic—yet plausible—definition. We measure the complexity of the model in terms of its entropic deviation from the simplest possible model: a network of isolated nodes. We know that the entropy of N isolated nodes is given by Siso = logN. This gives an upper bound on the von Neumann entropy of a hypergraph, guaranteeing that 0≤Sτ[d]≤Siso. Therefore, we define the model complexity C as:

C(ρτ[d]) = Siso − Sτ[d] ≥ 0: (23)

In fact, by definition we have that the complexity can take values in C(ρτ[d]) ∈ [0, Siso], reaching its minimum for hypergraphs with maximal entropy (e.g., isolated nodes or hypergraphs with high regularity), and its maximum for hypergraphs with minimal entropy (e.g., hypergraphs with more complex structures). From the point of view of the eigenvalues of the density matrix, the model complexity is expected to be lower when the eigenvalues of ρτ[d] are all similar, and larger when they are more diverse.

Minimizing the cost function
We can now define the cost function L by combining the information loss in Eq. (21) and the complexity in Eq. (23):

L(ρτ[dmax]|ρτ[d]) = DKL(ρτ[dmax]|ρτ[d]) + C(ρτ[d]) : (24)

By definition, minimizing the cost function corresponds to maximizing the accuracy of the model and, at the same time, minimizing the model complexity (Occam’s Razor). To obtain the best compression, we find the smallest order d which gives

dopt = min L(ρτ[dmax]|ρτ[d]) : (25)

It is worth mentioning that the propagation scale τ works as a resolution parameter. When τ is very large, the process approaches the steady state, and the network topology becomes irrelevant and, consequently, any model can be a good model. Whereas, at very small τ, the field evolution is linear and through the paths of length ≈ 1, exhibiting maximum resolution. Although we explore a variety of values of τ, we mainly focus on a characteristic propagation scale τc, the largest τ for which a linearization of the time evolution operator is still valid. Assume the eigenvalues of the Laplacian are given by {λℓ} where ℓ = 1, 2, . . . N and let λN be the largest of them. Then, the eigenvalues of the time-evolution operator e−τL are given by {e−τλℓ}. Here, τc = 1/λN, ensuring that the last eigenvalue of the time evolution operator is reasonably linearizable e−τcλN ≈ 1−τcλN. This ensures an acceptable linearization for the rest of the eigenvalues since λN is the largest eigenvalue of the Laplacian.

It is important to remark that in the original message length formulation the hypothesis is that observational data can be described by a generative model characterized by a probability distribution. In turn, such a distribution depends on some unknown parameters that can be fitted by maximizing the Bayesian posterior or, equivalently, minimizing the message length. Consequently, the term accounting for model complexity, C, depends on the number of parameters used to describe that probability distribution. Here, we deviate from the standard Bayesian approach, because we work with operators and not with probability distributions, thus a direct generalization of the Bayes theorem in this context is difficult. In fact, our model is a density matrix that characterizes the probability distribution of activating pathways for information flows60.

Consequently, using our formalism, there is no direct way to reconcile our definition of cost function with the standard Bayesian complexity term. In fact, the latter can be understood as Shannon’s surprise about the prior probability, which has no direct counterpart in our formalism. To partially overcome this issue, we use again a Kullback-Leibler divergence to approximate the model complexity and define an effective regularization term. In fact, this reduces to the deviation of Von Neumann entropy of the d-th order model from the case of a “gas network” (i.e., N disconnected nodes), which can be understood as the maximally mixed state to define our network. In fact, if we define ρiso = I/N, then DKL(ρ|ρiso) ≡ C encodes the expected “excess surprise” that we have about the state ρ if we use ρiso as a prior model for it:

DKL(ρ|ρiso) = Tr(ρ(logρ − logρiso)) = Tr(ρlogρ) − Tr(ρlogρiso) = −S − Tr(ρ(logI − log(N))) = logN − S ≡ C: (26)

Remarkably, the term logN is constant and does not alter the optimization procedure, but it is useful to reconnect—at least in spirit, and not formally—our complexity term C to the surprise about the prior appearing from the classical Bayesian approach. For this reason, throughout this paper, we will use the term cost function instead of description length to acknowledge the existing gap.

In this context, alternative definitions can be adopted for both the information loss and the complexity terms. We provide one such example in the Supplementary Information Sec. V and show that, when capturing different features, the resulting outcomes may differ.

Rescaling τ
Complete hypergraph. In some extremely regular structures, such as complete hypergraphs or some simplicial complex lattices, the Laplacians at all orders are proportional. For example, for complete hypergraphs,

L(d) = (K(d) / (N−1)) L(1), (27)

and consequently

L[d] = (d / (N−1)) L[1] = d L[1] : (28)

Another direct but useful consequence of this is the relationship between the multiorder Laplacians at any two orders

L[d] = (d/d0) L[d0] : (29)

Since by definition Eq. (3), the density matrix ρτ[d] depends only on the product τL[d], we can write

ρτ[d] = e−τdL[1] / Tr(e−τdL[1]) = ρeτ[1] (30)

where eτ=dτ, or equivalently, between two orders

ρτ[d] = ρ(d/d0)τ[d0] : (31)

Hence, instead of selecting a single diffusion time τ for all orders, we can select a different, more appropriate diffusion time at each order. Specifically, we can select a main τ, and rescale it at each order to ensure that the density matrices are the same. One could choose the rescaling so that the density matrices would be equal to any of them non-rescaled, e.g., to ρτ[1]. However, to ensure that the information loss vanishes when considering all possible orders, DKL(ρτ[dmax]|ρτ[dmax]) = 0, we need to set all of them equal to ρτ[1]. This is achieved by rescaling the diffusion time by

τ0(d) = (dmax / d) τ (32)

so that

ρτ0[d] = ρ(dmax/d)τ[d] = ρτ[dmax] : (33)

General case of proportional Laplacians. In general, the Laplacian of order d is proportional to that of order 1, L(d) ∝ L(1), when their respective adjacency matrices are proportional, that is

A(d) = d B(d) A(1) : (34)

where B(d) is a coefficient that may depend on the order d. Indeed, by definition, the generalized degree and adjacency matrix are related by K(d)i = (1/d) Σj A(d)ij, and thus we also have

K(d) = B(d)K(1), (35)

ensuring that the Laplacian matrices are proportional, L(d)=B(d)L(1). Equation (35) implies that

B(d) = 〈K(d)〉 / 〈K(1)〉, (36)

and hence, the multiorder Laplacian up to order d is given by

L[d] = (d / 〈K(1)〉) L(1) = d L[1] (37)

which is consistent with the complete hypergraph case in Eq. (27), where we have 〈K(1)〉 = N−1. The rest of the derivation is thus the same as in the complete graph case: rescaling τ per order as

τ0(d) = (dmax / d) τ (38)

ensures

ρτ0[d] = ρ(dmax/d)τ[d] = ρτ[dmax] : (39)

In turn, this yields no information loss and the flat cost function is given by L(ρτ[dmax]|ρτ[d0]) = C(ρτ[dmax]) ∀d.

Note again that, in general, a hypergraph does not need to have hyperedges at every order below D, unless it is a simplicial complex. If there is no hyperedge at some orders d, Eq. (37) must be adjusted, as the factor d comes from the number of orders present. The set of orders present is in general D={i≤d: 〈K(i)〉>0}, so that Eq. (37) becomes L[d]=|D|L[1] and Eq. (38) becomes τ0(d)=(dmax/|D|)τ.

Higher-order lattices. For a triangular lattice in which every triangle is promoted to a 2-simplex, each node is part of six 1-simplices and six 2-simplices. Furthermore, each 1-simplex of the lattice is part of two different 2-simplices. This means that each pair of nodes share two 2-simplices if they share one 1-simplex, and zero otherwise. Formally:

K(2)i = 6 = K(1)i and A(2)ij = 2A(1)ij, (40)

so that B(2)=1.

Description of the datasets
We assigned a category to each of the 60 empirical datasets. All datasets are accessible via XGI26 and stored at https://zenodo.org/communities/xgi/.

Coauthorship.—Each of the coauthorship datasets corresponds to papers published in a single year in geology (1980, 1981, 1982, 1983, 1984) and history (2010, 2011, 2012, 2013, 2014). A node represents an author, and a hyperedge represents a publication marked with the “Geology” or “History” tag in the Microsoft Academic Graph67.

Contact.—In the contact datasets, a node represents a person and a hyperedge represents a group or people in close proximity at a given time. Most of the original datasets are from the SocioPatterns collaboration68,69.

Drugs.—The drugs datasets include two constructed in ref. 24 with data from the National Drug Code Directory (NDC). In ndc-classes, a node represents a label (a short text description of a drug’s function) and a hyperedge represents a set of those labels applied to a given drug. In ndc-substances, a node represents a substance and a hyperedge is the set of substances in a given drug. The Drug Abuse Warning Network (DAWN) is a national health survey system that records drug use that contributes to hospital emergency department visits throughout the United States. A node represents a drug, and a hyperedge is the set of drugs used by a given patient (as reported by the patient) in an emergency department visit. For a period of time, the recording system only recorded the first 16 drugs reported by a patient, so the dataset only uses the first 16 drugs (at most).

Ecology.—In the ecology datasets, a node represents a plants species, and a hyperedge represents the set of pollinator species that visit a given plant species. These datasets are from the Web of Life (www.web-of-life.es).

Email.—In the email datasets, a node represents an email address and a hyperedge is the set of all recipient addresses included in an email, including the sender’s.

Metabolism.—Each metabolism dataset is a metabolic reaction networks at the genome scale, from the Biochemical, Genetic and Genomic (BiGG) models database (http://bigg.ucsd.edu/). A node represent a metabolite, and a hyperedge represents a set of metabolites participating in a metabolic reaction. Different datasets represents different metabolic reaction networks in different organisms. These datasets were analyzed as directed hypergraphs in ref. 70, but we converted them to unweighted undirected hypergraphs.

Politics.—In the politics datasets, a node represents a member of the US Congress/House/Senate and a hyperedge is the set of members co-sponsoring a bill between 1973 and 2016. In the House Committee dataset, a node is a member of the US House of Representatives, and a hyperedge is a set of member forming a committee (1993–2017). All datasets were constructed in ref. 24.

Tags.—In the tags datasets, a node represents a tag, and a hyperedge is a set of tags associated to a question on online Stack Exchange forums (Mathematics Stack Exchange and Ask Ubuntu). The tag datasets were constructed in ref. 24 with data from the Stack Exchange data dump.

Theater.—Each dataset is a Shakespeare theater play, where a node represents a character and a hyperedge represents a set of characters that appear together on stage in a scene. The Hyperbard datasets were introduced in ref. 71.

Randomizing strategies
We applied three randomization strategies: the configuration model, the node swap, and the shuffling.

The configuration model we used is a generalization of the standard pairwise configuration model to higher-order interactions25,59. We applied a configuration model to each order of the hypergraph separately. At each order, the configuration model connects stubs so as to randomize hyperedges while preserving the degree sequences at each order.

The node swap procedure was initially proposed in ref. 17. The idea is to preserve the degree distributions at each order, but to lower the cross-order degree correlations. To do this, at each order d>1, we pair nodes (i,j) and swap their generalized degrees K(d)i ↔ K(d)j by swapping their hyperedge memberships at that order. At each order, each node is paired with a single other nodes so that all nodes are swapped. To more efficiently reduce correlations between orders, high-degree nodes are paired with low-degree degree nodes in priority. Also, to avoid correlation between orders, the node pairing is randomly shuffled at each order.

The shuffling procedure is the most destructive one: it only preserves the number of hyperedges. It works by replacing each d-hyperedge by another, randomly selected, non-existing d-hyperedge. This procedure destroys intra- and inter-order correlations between hyperedges and tends to produce structures akin to random hypergraphs.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I first want to thank and congratulate the authors for addressing my concerns with such care. Since the original version, the\nauthors have made significant improvements to the manuscript and they have resolved my major concerns. In particular,\nthey have :\n- compiled the largest dataset of empirical higher-order networks in the field of Network Science and related fields (quality of\nthe data);\n- strengthened the narrative by reorganizing the text, adding additional remarks, refining the terminology, and including new\nreferences and figures;\n- improved the precision of their notation (e.g., in Eqs.~(5) and (6));\n- added a very insightful analytical result using the hyperring (lines 401 to 425 are particularly illuminating), significantly\nimproving the interpretation of the results in all the paper;\nThe results clarify how significant are higher-order interactions in complex systems and is definitely moving the field forward\n(see also my previous comments). After making the minor corrections below, especially meant to improve the flow of the\nwriting, I believe the work is perfectly suited for Nature Communications.\nMinor corrections\n————————\n1 - Below Eq. (6), d_opt^(\\tau) = d_max and d_opt^(\\tau) = 1.\n2 - Fig 3 caption : it is negative (positive) when dopt = 2 (1) -> dopt = 2 (= 1) to avoid having (1) like equation (1)\n3 - Fig 3(b), I did not see a description of the dashed lines and tau = 1 could simply be to the bottom right of the green curve\nto avoid cluttering\n4 - 329: \\forall -> for all\n5 - Eq. (11), ->\\left(dmax/d - 1 \\right) or \\Big( \\Big). Suggestion: it would be perhaps more easily readable if, after the equality,\nyou skip the line // so that all terms enter in the second line.\n6 - Lines 368 is hard to read. In fact, line 363 to 369 should be written more cleanly. It would probably help to write the\nexpectation E_rho… and lambda_k^delta side by side in one equation \\begin{align*} \\end{align*} after ‘where’ in line 363,\nfollowing with there respective description.\n7 - Lines 391-392 refer to Eqs. (11-12) specifically\n8 - Line 405: the second expectation should be delta instead of d\n9 - There is a contradiction between line 419 and Fig. 3(c) caption. In the former (which I believe is the correct one), it is said\nthat the hyperring is reducible at short time scales while in the latter, it is written that at tau_short (where the difference\nbetween cost function is negative), d_opt = 2 and it is irreducible. Suggestion: 3(c) could be made clearer by increasing its\nsize relative to 3(a,b) (that could be above or on the left side), drawing a grey zone where the difference is negative and\nwriting that chi_tau(H) = 0 in or above this zone.\n10 - Line 474: in simplicial complex, ‘where higher-order hyperedges are strongly nested within lower-order ones’ -> a bit\nconfusing, the sentence is heavy\n11 - Line 476: add little new spectral ...", "reviewer_2": "Many thanks again for the authors for this substantial round of revisions. I find the new framing of mechanisms vs behaviours\nmuch more suitable to the contribution of the paper, and the new analytical hyperring example very illuminating.\nI think this manuscript has improved significantly since the original submission, and now I am happy to support its\npublication in Nature Communications.\nHowever, I must say that two of my concerns still stand (albeit partially addressed):\n1) Information flow: I appreciate the authors' efforts to tone down some of their earlier claims around technical terms like\n\"message length\" and \"information flow\". I think this is now more representative of the work itself. Nonetheless, and despite\nthe author's justification around line 155, I think using \"information flow\" in the title is misleading as it does not involve\ninformation flow as used in information theory.\n2) Comparisons with Amari's work: I'm afraid I disagree with the somewhat restrictive view of the authors in their reply to my\n\"Comment 6\" (as numbered in the response document). Specifically, I disagree with the statement that \"these random variables can typically be seen as representing a multivariate time series.\" As an example of a common application of these\nmeasures, one could use a high-order Ising model (with a corresponding network structure), obtain its corresponding p(x),\nand calculate Amari's measure for each interaction order. This is directly comparable to the authors' work, but for discrete\nsystems. I don't think numerical experiments are necessary, but some discussion of this and other related measures coming\nfrom information theory is warranted."}, "author_response": "Response to Reviewer Reports\nNCOMMS-25-87411-T\nReducibility of higher-order networks via information flow\nWe are grateful to both Reviewers for their positive words and remaining minor comments.\nWe have revised the manuscript accordingly, and provide point-by-point responses below.\nResponse to Reviewer 1\nComment 1. I first want to thank and congratulate the authors for addressing\nmy concerns with such care. Since the original version, the authors have made sig-\nnificant improvements to the manuscript and they have resolved my major concerns.\nIn particular, they have :\n- compiled the largest dataset of empirical higher-order networks in the field of Net-\nwork Science and related fields (quality of the data);\n- strengthened the narrative by reorganizing the text, adding additional remarks, re-\nfining the terminology, and including new references and figures;\n- improved the precision of their notation (e.g., in Eqs. (5) and (6));\n- added a very insightful analytical result using the hyperring (lines 401 to 425 are\nparticularly illuminating), significantly improving the interpretation of the results in\nall the paper;\nThe results clarify how significant are higher-order interactions in complex systems\nand is definitely moving the field forward (see also my previous comments). After\nmaking the minor corrections below, especially meant to improve the flow of the\nwriting, I believe the work is perfectly suited for Nature Communications.\nReply:\nWethanktheReviewerfortheirfeedbackandpositivewords. Theirfeedbacksincethefirst\nround of reviews has pushed us and results in an improved manuscript. We have updated\nthe manuscript to address the minor corrections suggested.\nMinor corrections\n————————\n1 - Below Eq. (6), d_opt^(\\tau) = d_max and d_opt^(\\tau) = 1.\n2 - Fig 3 caption : it is negative (positive) when dopt = 2 (1) -> dopt = 2 (= 1)\nto avoid having (1) like equation (1)\n3 - Fig 3(b), I did not see a description of the dashed lines and tau = 1 could simply\nbe to the bottom right of the green curve to avoid cluttering\n4 - 329: \\forall -> for all\n5 - Eq. (11), ->\\left(dmax/d - 1 \\right) or \\Big( \\Big). Suggestion: it\nwould be perhaps more easily readable if, after the equality, you skip the line // so\nthat all terms enter in the second line.\n6 - Lines 368 is hard to read. In fact, line 363 to 369 should be written more cleanly.\nIt would probably help to write the expectation E_rho... and lambda_k^delta side\nby side in one equation \\begin{align*} \\end{align*} after ‘where’ in line 363,\nfollowing with there respective description.\n7 - Lines 391-392 refer to Eqs. (11-12) specifically\n8 - Line 405: the second expectation should be delta instead of d\n9 - There is a contradiction between line 419 and Fig. 3(c) caption. In the former\n(which I believe is the correct one), it is said that the hyperring is reducible at short\ntime scales while in the latter, it is written that at tau_short (where the difference\nbetween cost function is negative), d_opt = 2 and it is irreducible. Suggestion: 3(c)\ncould be made clearer by increasing its size relative to 3(a,b) (that could be above or\non the left side), drawing a grey zone where the difference is negative and writing\nthat chi_tau(H) = 0 in or above this zone.\n10 - Line 474: in simplicial complex, ‘where higher-order hyperedges are strongly\nnested within lower-order ones’ -> a bit confusing, the sentence is heavy\n11 - Line 476: add little new spectral ...\nReply:\nWe have fixed all of these, thank you for the careful reading. For 8., it was actually a typo\nand should have been d everywhere above—we fixed this. For 9., we fixed the text so that\neverything is consistent. We also clarified the Figure 3.\nResponse to Reviewer 2\nrevisions. I find the new framing of mechanisms vs behaviours much more suitable\nto the contribution of the paper, and the new analytical hyperring example very\nilluminating.\nI think this manuscript has improved significantly since the original submission, and\nnow I am happy to support its publication in Nature Communications.\nReply:\nWe thank the Reviewer for their constructive feedback since the first round of reviews and\nfor their kind words. We have addressed the remaining minor comments.\nComment 2. However, I must say that two of my concerns still stand (albeit\npartially addressed):\n1) Information flow: I appreciate the authors’ efforts to tone down some of their\nearlier claims around technical terms like ”message length” and ”information flow”.\nI think this is now more representative of the work itself. Nonetheless, and despite\nthe author’s justification around line 155, I think using ”information flow” in the\ntitle is misleading as it does not involve information flow as used in information\ntheory.\nReply:\nToavoidanyconfusion,wehavechangedthetitleto: “Reducibilityofhigher-ordernetworks\nfrom dynamics”.\nComment 3. 2) Comparisons with Amari’s work: I’m afraid I disagree with\nthe somewhat restrictive view of the authors in their reply to my ”Comment 6”\n(as numbered in the response document). Specifically, I disagree with the statement\nthat ”these random variables can typically be seen as representing a multivariate\ntime series.” As an example of a common application of these measures, one could\nuse a high-order Ising model (with a corresponding network structure), obtain its\ncorresponding p(x), and calculate Amari’s measure for each interaction order. This\nis directly comparable to the authors’ work, but for discrete systems. I don’t think\nnumerical experiments are necessary, but some discussion of this and other related\nmeasures coming from information theory is warranted.\nWe thank the Reviewer for this precision. We have now briefly mentioned this.\n4"}]

Decision letter