Ecological Network Structure is Multidimensional: Rethinking the Structure-Stability Debate

1 Introduction

Ecological theory has long sought general principles linking the architecture of interaction networks to ecosystem stability. Early work framed this relationship through system complexity, most prominently in Robert May’s analysis of random community matrices, which showed that increasing species richness, interaction density, and interaction-strength variance can reduce local dynamical stability in large ecological systems (May 1972, 1974). Subsequent work showed that particular forms of network organisation—including interaction structure, trophic hierarchy, and compartmentalisation—can alter stability outcomes (McCann 2000; Allesina and Tang 2012; Stouffer and Bascompte 2010; Tang et al. 2014). More recently, structural theories such as trophic coherence have emphasised that food webs occupy constrained regions of structural space in which particular architectural features influence stability (Johnson et al. 2014).

Despite these advances, there is little consensus on which aspects of network structure predict ecological stability. Different studies identify different structural predictors, sometimes with apparently conflicting effects. These inconsistencies are often interpreted as evidence that structure–stability relationships are weak, context-dependent, or system-specific. Here, we consider a complementary explanation - both network structure and stability are multidimensional, so different studies may be examining different parts of the same broader relationship.

Food webs are described using many structural metrics, including connectance, trophic organisation, modularity, omnivory, pathway structure, and motif composition. These metrics often covary, capture overlapping information, or describe organisation at different scales (Vermaat et al. 2009; Lau et al. 2017; Thompson et al. 2012). Rather than treating each metric as an independent candidate predictor, we therefore view them as observations of an underlying structural space. In this view, correlated metrics can be grouped into structural modules (coherent domains of network organisation) while multivariate axes such as principal components describe the dominant gradients through that space. A single metric is thus consequently a projection of structural space, not necessarily a universal descriptor of network structure.

Stability is similarly multidimensional. Ecological stability encompasses distinct responses such as resistance to disturbance, persistence or recovery following perturbation, and other properties of system dynamics (Domínguez-García et al. 2019; Ives and Carpenter 2007; Loreau and de Mazancourt 2013; Donohue et al. 2016; Chen et al. 2024). These properties arise from different processes and need not respond in the same way to changes in network structure. Thus, there is no a priori reason to expect one structural metric—or even one structural dimension—to predict all components of stability.

This leads to a simple conceptual framework. Ecological networks occupy a multidimensional structural space, stability can be represented by multiple stability components, and the structure–stability relationship is a mapping between these spaces. Under this framework, apparently conflicting results can arise when different structural features are related to different stability components. A predictor may therefore be positively associated with one component and negatively associated with another without the two relationships being contradictory.

Our empirical stability measures extend beyond a single dynamical property because empirical food webs are generally observed as interaction networks rather than as time series or fully parameterised dynamical systems. We therefore consider complementary network-level responses to perturbation such as robustness to species loss, recovery/persistence following sequential extinctions, spectral-radius-based stability potential, and structural controllability (Dunne et al. 2002; Jonsson et al. 2015; Staniczenko et al. 2013; Liu et al. 2011). These measures should not be interpreted as interchangeable definitions of stability. Rather, they represent distinct properties through which network structure may influence persistence, resistance, recovery, or control.

This conceptual framework is illustrated in Figure 1. Structural metrics are organised into a small number of underlying dimensions (Panel B), stability comprises multiple components (Panel A), and each stability component can depend on a different projection of structural space (Panel C). The resulting mapping is therefore not a single structure–stability curve, but a multidimensional set of relationships that may include trade-offs among stability components.

Figure 1: Conceptual framework linking multidimensional network structure (S) to multidimensional ecological stability (Ω). Structural variation is organised along a small number of axes (Panel B), while stability comprises multiple components (Panel A). Each stability component depends on a different projection of structural space (Panel C), generating process-dependent relationships and trade-offs among stability outcomes.

Here we test this framework using a compilation of empirical food webs characterised by a broad suite of structural descriptors. We first ask whether structural metrics form statistically robust modules based on their covariation across food webs. We then ask whether these modules align with the dominant multivariate axes of structural variation. Finally, we test whether different structural representations predict different components of stability, and whether individual structural features change direction across those components. We use these analyses to evaluate the central hypothesis that structure–stability relationships are best understood as a mapping between multidimensional structural and stability spaces rather than as universal relationships between individual metrics.

2 Materials & Methods

2.1 Data Compilation

We compiled food web data from Mangal (Poisot et al. 2016), Web of Life (Fortuna et al. 2014), and the canonical networks used by Vermaat et al. (2009), resulting in a total of XX networks. All networks were treated as binary and characterised using a suite of XX structural metrics (see Table 1).

Table 1: An informative caption about the different network properties. We use a combination of metrics from both the original Vermaat et al. (2009) paper as well as including those that have been identified by Thompson et al. (2012) and have been linked to emerging ecosystem properties such as stability

Label	Definition	Structural interpretation	Reference
Basal	Proportion of taxa with zero vulnerability (no consumers).	Quantifies the proportion of species representing basal energy inputs to the network.
Top	Proportion of taxa with zero generality (no resources).	Describes the relative prevalence of terminal consumers in the network.
Intermediate	Proportion of taxa with both consumers and resources.	Captures the proportion of species participating in both upward and downward energy transfer.
Richness (S)	Number of taxa (nodes) in the network.	Describes network size.
Links (L)	Total number of trophic interactions (edges).	Describes interaction density independent of network size.
Connectance	\(L/S^2\), where \(S\) is the number of species and \(L\) the number of links	Measures the proportion of realised interactions relative to all possible interactions.	Dunne et al. 2002
L/S	Mean number of links per species.	Captures average interaction density per taxon.
Cannibal	Proportion of taxa with self-loops.	Quantifies the prevalence of cannibalistic interactions.
Herbivore	Proportion of taxa feeding exclusively on basal species.	Describes the representation of primary consumers.
Trophic level (TL)	Prey-weighted trophic level averaged across taxa.	Captures the vertical organisation of energy transfer.	Williams and Martinez (2004)
MaxSim	Mean maximum trophic similarity of each taxon to all others.	Quantifies functional similarity based on shared predators and prey.	Yodzis and Winemiller (1999)
Centrality	Sum of differences between the maximum node centrality and all node centralities (Freeman’s centrality).	Captures how unevenly influence or connectivity is distributed among species.	Estrada and Bodin (2008)
ChLen	Mean length of all food chains from basal to top taxa.	Describes the average number of steps in energy-transfer pathways.
ChSD	Standard deviation of food chain length.	Captures variability in pathway lengths.
ChNum	Log-transformed number of distinct food chains.	Quantifies the multiplicity of alternative energy pathways.
Path	Mean shortest path length between all species pairs.	Describes the average distance between taxa within the network.
Diameter	Maximum shortest path length between any two taxa.	Captures the largest network distance between species.
Omnivory	Proportion of taxa feeding on resources at multiple trophic levels.	Describes vertical coupling of energy channels.	McCann (2000)
Loop	Proportion of taxa involved in trophic loops.	Quantifies the prevalence of cyclic interaction pathways.
Prey:Predator	Ratio of prey taxa (basal + intermediate) to predator taxa (intermediate + top).	Describes the overall shape of the trophic structure.
Clust	Mean clustering coefficient.	Measures the tendency for taxa sharing interaction partners to also interact with each other.	Watts and Strogatz (1998)
GenSD	Normalised standard deviation of generality.	Captures heterogeneity in the number of resources per taxon.	Williams and Martinez (2004)
VulSD	Normalised standard deviation of vulnerability.	Captures heterogeneity in the number of consumers per taxon.	Williams and Martinez (2004)
LinkSD	Normalised standard deviation of total links per taxon.	Quantifies variation in species connectivity.
Intervality	Degree to which taxa can be ordered along a single niche dimension.	Measures the extent of niche ordering in trophic interactions.	Stouffer et al. (2006)
S1 (Linear chain)	Frequency of three-node linear chains (A → B → C) with no additional links.	Captures the prevalence of simple, unbranched energy-transfer pathways.	Stouffer et al. (2007) Milo et al. (2002)
S2 (Omnivory)	Frequency of three-node motifs forming a feed-forward loop (A → B → C, A → C).	Describes vertical coupling of trophic levels within small subnetworks.	Stouffer et al. (2007) Milo et al. (2002)
S4 (Direct competition)	Frequency of motifs where two consumers share a single resource (A ← B → C).	Describes the occurrence of shared-resource structures among consumers.	Stouffer et al. (2007) Milo et al. (2002)
S5 (Apparent competition)	Frequency of motifs where one consumer feeds on two resources (A → B ← C).	Captures the prevalence of shared-predator structures among resources.	Stouffer et al. (2007) Milo et al. (2002)
trophicVar	Measure of how much the trophic positions of species deviate from the mean trophic level.	Variance is linked to chain length, with low variance indicating few long chains and high variance indicates long chains and varied interactions	Pimm (1982)
TrophicCoherence	Measure how well the species in a food web fit into discrete trophic levels.	Highly coherent, neatly layered food webs are more stable to perturbations	Johnson et al. (2014)

In addition to network structural descriptors, we quantified four complementary stability-related properties of the food webs. First, robustness to species loss (\(R_{50}\)) quantifies the proportion of primary extinctions required to reduce network richness to 50% of its initial value, providing a threshold-based measure of tolerance to species removal (Dunne et al. 2002; Jonsson et al. 2015). Second, we quantified recovery/persistence using the area under the extinction curve, which relates primary species loss to secondary extinctions and captures how rapidly structural collapse propagates as perturbations accumulate. Third, spectral radius (\(\rho\)), defined as the largest real part of the eigenvalues of the adjacency matrix (Staniczenko et al. 2013), was used as a spectral measure of stability potential (Bascompte et al. 2003). Finally, structural controllability quantifies the minimum number of driver species required to control the system, capturing a distinct property of network reorganisation and directed influence (Liu et al. 2011).

We refer to these four outcomes as stability components for brevity, while recognising that they capture different processes rather than a single interchangeable quantity. They comprise (i) resistance to species loss, (ii) recovery/persistence under sequential species removal, (iii) spectral stability potential, and (iv) controllability. Because robustness and recovery/persistence depend on stochastic extinction sequences, we averaged each measure across 100 extinction runs.

Finally, we incorporated two complementary measures that capture higher-order aspects of network organisation. First, we quantified complexity based on the singular value spectrum of the adjacency matrix. Specifically, we computed the entropy of the singular values obtained via singular value decomposition (SVD), which captures how structural variance is distributed across orthogonal modes. Networks with more evenly distributed singular values exhibit higher complexity, reflecting greater heterogeneity in interaction pathways and reduced dominance of any single structural mode (Strydom et al. 2021). This provides a spectral analogue to traditional structural descriptors, linking network organisation to the distribution of interaction strength across latent dimensions.

Second, we included a composite scaling term of the form \(-\frac{1}{2}(\log{(Richness)} + \log{(Connectance)})\) This term captures the joint scaling of network size and density and can be interpreted as a log-transformed geometric mean of these quantities. This formulation is motivated by classical stability results, where system behaviour depends on the combined scaling of species richness and interaction density, rather than either quantity in isolation. As such, this term provides a compact representation of size–complexity tradeoffs and complements both local (e.g., degree distributions) and global (e.g., spectral) descriptors of network structure.

2.2 Identification of Structural Modules

We tested the hypothesis that structural descriptors of food webs are organised into statistically coherent modules reflecting shared ecological function or scaling relationships, rather than forming arbitrary clusters driven by sampling noise. If such modular organisation exists, then metrics within a module should exhibit strong internal correlation relative to metrics assigned to different modules. The resulting clusters should be robust to resampling of the data. The identified modular structure should exceed expectations under null models of random association among metrics.

We quantified pairwise associations among the XX structural metrics using Pearson correlations computed across food webs. Because ecological metrics may covary either positively or negatively depending on scaling relationships, we constructed two alternative distance matrices: \(1−r\), which preserves the sign of correlations and distinguishes positive from negative association and \(1−∣r∣\), which groups variables based on the magnitude of their association regardless of sign. These two definitions allow us to test whether modular structure depends on directional relationships or simply on the strength of coupling among metrics. Hierarchical clustering was performed using average linkage on each distance matrix to identify candidate structural modules.

The optimal number of clusters was evaluated across a range of partition sizes (\(k = 2–10\)) using average silhouette width to assess within-cluster cohesion and between-cluster separation. To further evaluate cluster robustness, we implemented bootstrap resampling with 1,000 bootstrap replicates. This procedure estimates approximately unbiased (AU) p-values for each cluster, quantifying the probability that a cluster is supported under repeated resampling of the data. Clusters were considered statistically robust when they exhibited high silhouette support, and approximately unbiased bootstrap support ≥ 0.95. This dual criterion ensures that identified modules are both structurally coherent and stable to sampling variation.

2.3 Multivariate Structure and Module–Axis Alignment

To evaluate whether structural metrics of food webs organise into coherent multivariate modules and whether those modules define principal axes of variation, we conducted a principal component analysis (PCA) followed by a permutation-based test of module–axis alignment. Under this framework, principal components define dominant structural gradients in the metric space, whereas modules represent hypothesised mechanistic groupings of structurally related metrics. Demonstrating alignment between these two representations would suggest that modular decomposition captures fundamental axes of ecological variation.

Skewed count-based metrics (e.g., link counts, interval counts, and related quantities) were log-transformed using \(log(x + 1)\) to reduce right skew. All remaining variables were standardised to zero mean and unit variance prior to analysis to ensure that metrics measured on different scales contributed equally to the ordination.

To quantify how structural modules align with principal axes, we decomposed variance in each principal component according to module membership. For each module \(m\) and principal component \(k\), we computed:

\[ A_{mk} = \sum_{i \in m} l^{2}_{ik} \]

where \(l_{ik}\) is the loading of metric \(i\) on \(PC_k\) and \(A_{mk}\) represents the fraction of variance in \(PC_k\) attributable to metrics in module \(m\). Because squared loadings sum to one within each principal component, this provides a direct partitioning of PC variance across modules. This produces a module × PC matrix describing geometric alignment between modular structure and multivariate axes.

To evaluate whether observed module–PC alignment exceeded expectations under random module structure, we implemented a permutation test. We randomly permuted module labels among metrics 1 000 times while preserving the number and size distribution of modules, and the PCA loadings. For each permutation \(r\), we recomputed \(A^{(r)}_{mk}\) so as to form a null distribution for each module–PC pair. From this we computed p-values and corresponding z-scores, which was then used to infer significant alignment when \(p_{mk} < 0.05 \land |Z_{mk}| > 1.96\)

where:

\[ p_{mk} = P(A^{(r)}_{mk} > A_{mk}) \]

and:

\[ Z_{mk} = \frac{A_{mk} - \mu_{mk}}{\sigma_{mk}} \]

Additionally we also evaluated if modules collectively aligned with PCA structure beyond random expectation. Overall concentration of variance within module–PC space was calculated as \(T = \sum_{m, k} A^{2}_{mk}\) and significance was determined by comparing the observed (\(T\)) the permutation distribution (\(T^{(r)}\)). Where \(p_{global} = P(T^{(r)} \geq T)\). A significant result indicates that module structure explains multivariate variance better than random groupings.

2.4 Representing the Multidimensional Structural Space

To evaluate how network structure predicts ecological stability, we first reduced a high-dimensional set of topological metrics into alternative, biologically interpretable predictor sets. Because dimensionality reduction can fundamentally shape inference, we explicitly compared two conceptually distinct representations of network structure: (1) structural domain representatives and (2) latent axes of network variation.

2.4.1 Structural-module representatives

We previously grouped network metrics into clusters based on their pairwise correlations, such that each cluster represented a structurally coherent domain (e.g., trophic composition, centralisation, path structure). Clustering was performed using correlation-based distance, ensuring that metrics grouped together reflected shared structural information rather than raw scale. To construct a reduced predictor set from these modules, we selected a single representative metric per module using a medoid approach. Within each cluster, we computed the absolute correlation matrix among member metrics and defined distance as \(1−∣r∣\). The medoid was identified as the metric minimising mean pairwise distance to other metrics in the cluster. This approach preserves structural diversity while minimising redundancy and does not rely on principal component loadings, which can bias selection toward dominant variance axes. The resulting ‘cluster medoid’ set retained one interpretable metric per structural module. Importantly, this procedure prioritises preservation of structural domains rather than overall variance magnitude, allowing low-variance but potentially mechanistically relevant features of network topology to be retained.

2.4.2 Latent structural axes

As a complementary representation of network structure, we performed principal component analysis (PCA) on the scaled full metric matrix. Principal components were retained until cumulative explained variance exceeded 80%, yielding a set of orthogonal axes describing the dominant gradients of variation in network topology. These retained PC scores constitute a low-dimensional latent representation of network space. Unlike cluster medoids, which preserve structural categories, PC scores preserve maximal variance and ensure orthogonality among predictors. We also explored a third reduction approach—selecting the single highest-loading metric per retained PC. However, because principal components are linear combinations of multiple metrics, this PC-dominant metric strategy substantially reduces the variance represented by each axis. Consequently, this approach captures less of the total structural variance.

2.4.3 Contrasting structural representations

The cluster-medoid and PCA-score representations reflect different theoretical assumptions about how structure relates to stability. The cluster-medoid approach assumes that ecological stability responds to discrete structural domains that may vary independently and need not align with the dominant gradients of network variation. In contrast, the PCA-score approach assumes that stability responds primarily to the major axes of variation in network topology, regardless of their interpretability in terms of structural domains. By explicitly comparing predictive performance across these representations, we test whether ecological stability is better explained by domain-specific structural mechanisms or by global geometric gradients of network organisation.

2.5 Linking Structural Space to Stability Components

To evaluate how hierarchical representations of network structure explain variation in stability, we modelled four complementary stability outcomes: Robustness (\(R_{50}\)), Resilience, Spectral Radius (\(\rho\)), and Structural Controllability. These outcomes represent distinct components of network response and should not be interpreted as interchangeable measures of a single stability property. Elastic net regression models were fit separately for each stability metric using three structural representations (cluster medoids, PC-dominant metrics, and full PCA scores) and XXX univariate, canonical proxies of stability (complexity, connectance, trophic coherence), allowing us to compare the structural determinants of different stability processes. All predictors and response variables were standardised prior to analysis to allow direct comparison of coefficient magnitudes.

Elastic net regularisation was chosen because it permits inference under correlated predictors by interpolating between ridge (distributed shrinkage) and lasso (sparse selection) penalties. The mixing parameter \(\alpha\) (0–1) determines the degree of sparsity, with lower \(\alpha\) values indicating distributed influence across many correlated structural descriptors (ridge-like behaviour), whereas higher \(\alpha\) values indicate that stability is driven by a smaller subset of structural predictors (lasso-like sparsity).

2.5.1 Model Evaluation

Predictive performance was assessed using repeated v-fold cross-validation (5 folds × 10 repeats) this improves the robustness of hyperparameter selection, particularly under correlated predictors and variable sample sizes. Repeating the partitioning reduces sensitivity to any single fold split and stabilises the selection of both \(\alpha\) and \(\lambda\), ensuring that the reported coefficients and variance contributions reflect consistent structural effects rather than idiosyncrasies of a particular data split. This approach enhances confidence that observed switching patterns in predictor importance across stability metrics represent genuine structural relationships rather than artefacts of sampling variability. Within each training partition, models were fit across a grid of \(\alpha\) values (0–1 in increments of 0.25), and for each \(\alpha\), the optimal penalty strength (\(\lambda\)) was selected via internal cross-validation on the training data. Specifically, \(\lambda\) values were chosen from a dense path generated by the elastic net algorithm, and the \(\lambda\) that maximised cross-validated \(R^2\) within the training fold was recorded. This process ensured that \(\lambda\) selection was independent of the held-out test data. Model performance was quantified as the mean cross-validated \(R^2\) computed on held-out test partitions across all repeats and folds.

After cross-validation, final models were refit to the full dataset using the mean selected \(\alpha\), and the λ value closest to the cross-validated optimum was used to extract standardised regression coefficients. These coefficients characterise the direction and relative magnitude of structural effects on each stability metric.

The primary objective was to characterise how structural organisation relates to different stability components rather than to identify a single optimal predictor set. The elastic net mixing parameter \(\alpha\) provides insight into the architecture of structural control, indicating whether stability is influenced by many predictors (ridge-like) or a sparse subset (lasso-like). To quantify the relative contribution of different structural modules, standardised regression coefficients from the final models were squared and aggregated within modules. The proportion of variance explained by each module was then scaled by the model’s cross-validated \(R^2\), yielding a measure of the absolute variance in stability explained by each structural component.

3 Results

3.1 Structural Metrics Organise into Robust Modules

Hierarchical clustering of structural descriptors revealed clear modular organisation within food web architecture Figure 2. Silhouette analysis indicated an optimal partition of the signed correlation matrix at k=7 modules, with a maximum average silhouette width of [INSERT VALUE]. Clustering based on absolute correlations produced a slightly coarser but comparable solution (k=5), demonstrating that the modular structure is robust to the treatment of correlation sign and reflects strong underlying association patterns among metrics.

Figure 2: Hierarchical clustering of food web structural descriptors. Dendrogram illustrating the modular organisation of structural metrics based on a signed correlation matrix. Clustering was performed using UPGMA. Coloured bars indicate the seven distinct modules, reflecting robust association patterns among metrics. The horizontal axis represents the correlation distance (1−r), where shorter branches indicate stronger associations between descriptors.

Bootstrap resampling further supported the stability of the major clusters. Most principal modules exhibited high approximately unbiased (AU) support values (AU ≥ [INSERT VALUE]), indicating that the identified clusters are not artefacts of sampling variability but represent statistically robust groupings of structural descriptors.

The resulting modular partition identified coherent groups of metrics describing distinct aspects of network architecture Figure 3. These included a large-scale architectural module integrating size and motif-based descriptors, modules capturing trophic organisation and energy routing, and more isolated structural descriptors that vary relatively independently from broader groupings. The seven-module solution Figure 3 reveals a hierarchical decomposition of food web structural space into interpretable ecological dimensions. Importantly, these modules can be understood not only as statistical groupings, but as reflecting successive constraints operating during network assembly. Comparison with both connectance and degree-constrained null models reinforces this interpretation by showing that not all structural dimensions arise from the same underlying processes (SUPP MATT REF). Early-emerging dimensions associated with network size and interaction density define the combinatorial space of possible interactions, consistent with classical complexity arguments in which system behaviour is governed by its size (May 1972; Allesina and Tang 2012). In contrast, higher-order modules capture constraints imposed by trophic organisation, energy flow, and species roles, reflecting ecological structuring processes described by niche-based and trophic theories of food web assembly (Williams and Martinez 2000; Stouffer and Bascompte 2010; Johnson et al. 2014). Collectively, this structure indicates that food web organisation arises from multiple, partially independent axes of variation that correspond to distinct constraints operating during network assembly.

Figure 3: Ecological network structure emerges through hierarchical constraints on assembly. Food web architecture can be decomposed into six structural modules organised across three levels of constraint. Combinatorial scaling defines the space of possible interactions as a function of system size. This baseline gives rise to statistical organisation, where degree constraints generate heterogeneity in species roles and regulatory structure (trophic role asymmetry, control heterogeneity). At the highest level, ecological structuring imposes trophic ordering, energy flow constraints, and the emergence of resource and control hubs (trophic integration, energy transport, flow control).

3.2 Structural Modules Align with Dominant Axes of Network Variation

Principal component analysis of standardised structural metrics revealed strong dimensional compression in food web structure Figure 4. The first three principal components accounted for 68.5% of total variance (PC1 = 31.2%, PC2 = 21.5%, PC3 = 15.6%), with 80.7% explained by the first five components, indicating that a small number of dominant gradients capture most structural variation.

Figure 4: PCA loadings plot illustrating the dominant gradients of food web variation. Each point denotes a structural metric, with colours corresponding to the modular groupings. The clustering of like-coloured points along specific PC axes demonstrates that the modular organisation is preserved within a reduced dimensional space.

Variance decomposition of principal components by module revealed a structured alignment between the modular organisation of metrics and the dominant axes of variation Figure 5. The first two principal components were jointly structured by Macro Complexity and Trophic Integration, but in contrasting ways indicating that variation primarily separates highly integrated trophic structure from size-driven architecture. PC3 represented a distinct structural dimension dominated by Control Heterogeneity, with additional contributions from Trophic Role Asymmetry, reflecting variation in the distribution of predation pressure and imbalances among trophic roles. This axis captures differences in how strongly interactions are concentrated among species and how unevenly trophic roles are distributed across the network.

Figure 5: Heatmap showing the correspondence between the seven identified food web modules (rows) and the principal components (columns) capturing 80% of cumulative variance. High values (darker shading) indicate strong alignment, where a specific module’s internal structure drives the variation captured by that PC.

Beyond the first three axes, additional modules aligned with more specialised components of structural variation. In particular, Transport Efficiency exhibited strong alignment with higher-order components (notably PC6), indicating that the geometry of trophic pathways varies largely independently of both network size and trophic integration. The Resource and Control Hubs module showed weaker and more diffuse alignment across higher-order components, suggesting that variation in energy entry points and centrality represents a secondary structural dimension.

A global permutation test confirmed that the observed alignment between modules and principal components was significantly stronger than expected under random assignment (\(p = 0.006\)), demonstrating that the modular structure captures meaningful geometric organisation in the multivariate space of network metrics Figure 6. Alignment was concentrated in a subset of principal components, indicating that modules map onto specific structural gradients rather than uniformly spanning the full space.

Figure 6: Heatmap showing module–principal component alignment quantified by permutation-derived Z-scores. Colour intensity reflects the magnitude and direction of alignment, and asterisks denote cells that meet joint significance criteria (permutation P < 0.05 and |Z| > 1.96). Rows correspond to structural metric modules and columns to principal components. Clustering highlights patterns of shared alignment across modules and multivariate structural gradients.

Collectively, these results demonstrate that food web structure is organised along a small number of dominant, ecologically interpretable axes. A primary structural gradient reflects a strong trophic organisation signal, while secondary axes capture variation in regulatory heterogeneity and structural scaling. This hierarchical organisation indicates that food web architecture is governed by multiple, partially independent dimensions of structural variation.

3.3 Mapping Structural Modules onto Stability Components

The alternative structural representations differed in variance retention, redundancy, and interpretability. Because these properties are methodological rather than ecological results, detailed comparisons among representations are provided in the Supplementary Material (Fig. S1). Here we focus on the resulting structure–stability relationships. The main comparison asks whether stability components are better understood through structural modules, dominant multivariate axes, or individual canonical predictors.

3.3.1 Different structural modules contribute to different stability components

3.3.1.1 Structural Modules Differentially Contribute to Stability Processes

Variance partitioning revealed that different stability components draw on different subsets of structural modules (Figure 7). Rather than identifying a single structural regime that consistently explains stability, the models showed component-specific patterns of structural contribution. Modules associated with trophic heterogeneity and role asymmetry contributed disproportionately to explaining Controllability and Recovery/Persistence, with metrics such as GenSD alone accounting for a large fraction of explained variance (up to ~50% in some models). In contrast, Stability Potential was primarily explained by modules capturing global network organisation, consistent with its dependence on spectral properties of the adjacency matrix. Resistance, however, showed relatively weak and diffuse contributions across modules, with no single domain exerting strong control. This reinforces the view that resistance is an emergent property arising from distributed features of network organisation, rather than a specific structural mechanism. Together, these patterns indicate that different stability components are associated with different regions of structural space.

Figure 7: Contribution of structural modules to ecosystem stability. Bars show the absolute variance in each stability metric explained by different structural modules, as estimated via elastic net regression. Facets correspond to the four stability metrics: Resistance (robustness, \(R_{50}\)), Recovery/Persistence (resilience AUC), Stability Potential (Spectral Radius (ρ)), and Controllability (Structural Controllability). Variance contributions are scaled by the cross-validated \(R^2\) of the corresponding model, allowing comparison across metrics. Differences in bar heights across facets indicate which modules are more influential for particular stability processes.

3.3.1.2 Structural effects reverse across stability components

Predictor effects frequently changed sign across stability components (Figure 8). Structural features associated positively with one component were often associated negatively with another, indicating that their effects depend on the stability property being considered. The most pronounced example is network complexity, which shows strong negative associations with both Resistance and Stability Potential, but becomes positively associated with Controllability. This indicates that increasing structural heterogeneity simultaneously reduces tolerance to perturbation while enhancing the system’s capacity for external control.

Similar reversals are observed for other structural features. For example, trophic composition variables such as the proportion of top predators consistently reduce resistance and stability potential, yet exhibit weaker or opposing effects in other stability contexts. Likewise, pathway-based metrics such as chain length strongly destabilise recovery dynamics while contributing differently to other components. Thus, structural features were not uniformly stabilising or destabilising across the four outcomes; instead, their associations were component-specific and could imply trade-offs among stability properties.

Figure 8: Standardised elastic net regression coefficients linking structural predictors to ecosystem stability metrics: Resistance (robustness, \(R_{50}\)), Recovery/Persistence (resilience AUC), Stability Potential (Spectral Radius (ρ)), and Controllability (Structural Controllability). Each point represents a predictor, with vertical lines extending from zero to the coefficient value. Coefficients are coloured by structural module and faceted by both network representation (Cluster Dominant, PC Dominant, PCA Axis, or Univariate) and stability metric. Positive coefficients indicate that higher values of the predictor increase the stability metric, while negative coefficients indicate a suppressive effect. Variation in coefficients across facets reveals switching of structural predictors’ influence on different stability outcomes.

3.3.1.3 Classical complexity captures only part of the mapping

Finally, comparison with classical univariate predictors highlights the limitations of traditional approaches. Metrics derived from complexity theory, including the May scaling term, explain only a small fraction of variation across stability components (typically \(R^2 < 0.1\)), and fail to capture the dominant structural drivers identified in multivariate models. While complexity shows strong negative associations with some stability components, consistent with classical predictions, it does not generalise across stability processes and cannot account for observed sign reversals. This indicates that classical complexity–stability relationships represent partial projections of a higher-dimensional system, capturing only one aspect of the broader mapping between network structure and ecological stability.

Figure 9: A vibe but don’t know where it belongs - supp mat and it beautifully mirrors the conceptual figure in the intro - just saying.

4 Discussion

4.1 A multidimensional structure–stability relationship

Our results support a multidimensional view of the relationship between food web structure and stability. Structural metrics did not behave as a collection of independent descriptors. Instead, they formed coherent modules that captured distinct aspects of network organisation, and these modules aligned with a small number of dominant multivariate axes. At the same time, the four stability components showed different structural associations. The resulting picture is therefore not a single relationship between ‘complexity’ and ‘stability’, but a mapping between multiple structural dimensions and multiple stability components.

This distinction provides a useful way to interpret apparently conflicting structure–stability relationships. A structural feature does not have to be intrinsically stabilising or destabilising: its association can depend on which component of stability is considered. We observed this directly in the frequent reversal of predictor effects across stability components. The same structural feature could be associated negatively with resistance or stability potential while being positively associated with controllability. Such reversals are difficult to reconcile within a one-dimensional framework, but are expected when different stability components depend on different projections of structural space.

4.2 Structural modules reveal organised dimensions of food web architecture

The clustering results indicate that commonly used network metrics can be reduced into a smaller set of coherent structural modules. These modules should be distinguished from the principal components used to describe the geometry of the same metric space. Modules group metrics according to their covariation and provide interpretable structural domains; principal components instead identify orthogonal axes that capture dominant variation. Their alignment in our data indicates that the modular organisation is not simply an arbitrary clustering of correlated variables, but is represented in the dominant gradients of network structure.

The null-model comparisons further suggest that not all structural organisation can be explained by generic link-level constraints. Modules associated with network size and degree structure were more readily reproduced under constrained randomisation, whereas other modules showed weaker correspondence with null expectations. This is consistent with the possibility that different dimensions of food web architecture arise from different constraints during network assembly. We treat this as an interpretation of the structural organisation rather than a direct test of assembly mechanisms; identifying the processes that generate particular modules remains an important direction for future work.

The classical complexity perspective provides a useful special case within this broader framework. The May scaling term showed negative associations with some stability components, consistent with the direction of classical complexity–stability expectations, but explained relatively little variation across the full set of outcomes. We therefore do not interpret these results as a rejection of complexity–stability theory. Instead, they suggest that classical complexity captures one projection of a broader structure–stability relationship and cannot, by itself, represent the full set of structural dimensions relevant to different stability components.

4.3 Implications for ecological theory and network modelling

A multidimensional view has implications for how ecological networks are compared, reconstructed, and modelled. If empirical networks occupy a structured but relatively low-dimensional region of structural space, then the ability to reproduce a single global property such as connectance is unlikely to guarantee that a synthetic network captures the structural features relevant to a particular stability component. Conversely, preserving multiple structural modules may provide a more informative target for network reconstruction and model evaluation.

The results also suggest that the choice of structural metrics should be driven by the ecological question rather than by the search for a universal set of predictors. Metrics that belong to the same structural module may provide largely overlapping information, whereas metrics from different modules may capture complementary dimensions. A principled strategy is therefore to select descriptors that span the structural modules relevant to the stability component or ecological process of interest while limiting unnecessary redundancy.

Finally, the component-specific associations highlight the importance of treating stability as a family of related properties rather than a single outcome. Apparent trade-offs among stability components may arise because the same structural organisation has different consequences for resistance, recovery/persistence, stability potential, and controllability. Understanding these trade-offs will require linking structural dimensions to explicit ecological and dynamical mechanisms, rather than assuming that a single metric can stand in for stability as a whole.

5 Conclusion

Food web structure is multidimensional, but much of its variation can be described by a relatively small number of coherent structural dimensions. Stability is likewise multidimensional, with different components associated with different regions of structural space. Our results therefore support a framework in which structure–stability relationships are understood as mappings between these two spaces rather than as universal effects of individual network metrics.

Within this framework, apparently conflicting relationships can arise because different structural features influence different stability components, and because the same feature can have opposing associations across components. Classical complexity–stability relationships remain informative as specific projections of this broader system, but they do not capture the full multidimensional mapping. A more complete understanding of ecological stability will therefore require approaches that explicitly represent structural dimensions, distinguish among stability components, and test the trade-offs that emerge between them.

References

Allesina, Stefano, and Si Tang. 2012. “Stability Criteria for Complex Ecosystems.” Nature 483 (7388, 7388): 205–8. https://doi.org/10.1038/nature10832.

Bascompte, J., P. Jordano, C. J. Melian, and J. M. Olesen. 2003. “The Nested Assembly of Plant-Animal Mutualistic Networks.” Proceedings of the National Academy of Sciences 100 (16): 9383–87. https://doi.org/10.1073/pnas.1633576100.

Chen, Can, Xu-Wen Wang, and Yang-Yu Liu. 2024. “Stability of Ecological Systems: A Theoretical Review.” Physics Reports, Stability of Ecological Systems: A Theoretical Review, vol. 1088 (October): 1–41. https://doi.org/10.1016/j.physrep.2024.08.001.

Domínguez-García, Virginia, Vasilis Dakos, and Sonia Kéfi. 2019. “Unveiling Dimensions of Stability in Complex Ecological Networks.” Proceedings of the National Academy of Sciences 116 (51): 25714–20. https://doi.org/10.1073/pnas.1904470116.

Donohue, Ian, Helmut Hillebrand, José M. Montoya, et al. 2016. “Navigating the Complexity of Ecological Stability.” Ecology Letters 19 (9): 1172–85. https://doi.org/10.1111/ele.12648.

Dunne, Jennifer, Richard J Williams, and Néo D Martinez. 2002. “Network Structure and Biodiversity Loss in Food Webs: Robustness Increases with Connectance.” Ecol. Lett. 5 (4): 558–67.

Estrada, Ernesto, and Örjan Bodin. 2008. “Using Network Centrality Measures to Manage Landscape Connectivity.” Ecological Applications 18 (7): 1810–25. https://doi.org/10.1890/07-1419.1.

Fortuna, Miguel A., Raul Ortega, and Jordi Bascompte. 2014. “The Web of Life.” Pre-published March 11. https://doi.org/10.48550/arXiv.1403.2575.

Ives, Anthony R., and Stephen R. Carpenter. 2007. “Stability and Diversity of Ecosystems.” Science 317 (5834): 58–62. https://doi.org/10.1126/science.1133258.

Johnson, Samuel, Virginia Domínguez-García, Luca Donetti, and Miguel A. Muñoz. 2014. “Trophic Coherence Determines Food-Web Stability.” Proceedings of the National Academy of Sciences 111 (50): 17923–28. https://doi.org/10.1073/pnas.1409077111.

Jonsson, Tomas, Sofia Berg, Alexander Pimenov, Catherine Palmer, and Mark Emmerson. 2015. “The Reliability of R50 as a Measure of Vulnerability of Food Webs to Sequential Species Deletions.” Oikos 124 (4): 446–57. https://doi.org/10.1111/oik.01588.

Lau, Matthew K., Stuart R. Borrett, Benjamin Baiser, Nicholas J. Gotelli, and Aaron M. Ellison. 2017. “Ecological Network Metrics: Opportunities for Synthesis.” Ecosphere 8 (8): e01900. https://doi.org/10.1002/ecs2.1900.

Liu, Yang-Yu, Jean-Jacques Slotine, and Albert-László Barabási. 2011. “Controllability of Complex Networks.” Nature 473 (7346): 167–73. https://doi.org/10.1038/nature10011.

Loreau, Michel, and Claire de Mazancourt. 2013. “Biodiversity and Ecosystem Stability: A Synthesis of Underlying Mechanisms.” Ecology Letters 16 (s1): 106–15. https://doi.org/10.1111/ele.12073.

May, Robert M. 1972. “Will a Large Complex System Be Stable?” Nature 238 (5364): 413–14. https://doi.org/10.1038/238413a0.

May, Robert M. 1974. Stability and Complexity in Model Ecosystems. Vol. 1. Princeton University Press. https://doi.org/10.2307/j.ctvs32rq4.

McCann, Kevin Shear. 2000. “The Diversity–Stability Debate.” Nature 405 (6783): 228–33. https://doi.org/10.1038/35012234.

Milo, R., S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon. 2002. “Network Motifs: Simple Building Blocks of Complex Networks.” Science 298 (5594): 824–27. https://doi.org/10.1126/science.298.5594.824.

Pimm, Stuart L. 1982. Food Webs. Springer Netherlands. https://doi.org/10.1007/978-94-009-5925-5.

Poisot, Timothée, Benjamin Baiser, Jennifer Dunne, et al. 2016. “Mangal – Making Ecological Network Analysis Simple.” Ecography 39 (4): 384–90. https://doi.org/10.1111/ecog.00976.

Staniczenko, Phillip P. A., Jason C. Kopp, and Stefano Allesina. 2013. “The Ghost of Nestedness in Ecological Networks.” Nature Communications 4 (1): 1391. https://doi.org/10.1038/ncomms2422.

Stouffer, Daniel B., and Jordi Bascompte. 2010. “Understanding Food-Web Persistence from Local to Global Scales.” Ecology Letters 13 (2): 154–61. https://doi.org/10.1111/j.1461-0248.2009.01407.x.

Stouffer, Daniel B., Juan Camacho, and Luís A. Nunes Amaral. 2006. “A Robust Measure of Food Web Intervality.” Proceedings of the National Academy of Sciences 103 (50): 19015–20. https://doi.org/10.1073/pnas.0603844103.

Stouffer, Daniel B, Juan Camacho, Wenxin Jiang, and Luís A Nunes Amaral. 2007. “Evidence for the Existence of a Robust Pattern of Prey Selection in Food Webs.” Proceedings of the Royal Society B: Biological Sciences 274 (1621): 1931–40. https://doi.org/10.1098/rspb.2007.0571.

Strydom, Tanya, Giulio V. Dalla Riva, and Timothée Poisot. 2021. “SVD Entropy Reveals the High Complexity of Ecological Networks.” Frontiers in Ecology and Evolution 9. https://doi.org/10.3389/fevo.2021.623141.

Tang, Si, Samraat Pawar, and Stefano Allesina. 2014. “Correlation Between Interaction Strengths Drives Stability in Large Ecological Networks.” Ecology Letters 17 (9): 1094–100. https://doi.org/10.1111/ele.12312.

Thompson, Ross M., Ulrich Brose, Jennifer Dunne, et al. 2012. “Food Webs: Reconciling the Structure and Function of Biodiversity.” Trends in Ecology & Evolution 27 (12): 689–97. https://doi.org/10.1016/j.tree.2012.08.005.

Vermaat, Jan E., Jennifer A. Dunne, and Alison J. Gilbert. 2009. “Major Dimensions in Food-Web Structure Properties.” Ecology 90 (1): 278–82. https://doi.org/10.1890/07-0978.1.

Watts, Duncan J., and Steven H. Strogatz. 1998. “Collective Dynamics of ‘Small-World’ Networks.” Nature 393 (6684): 440–42. https://doi.org/10.1038/30918.

Williams, Richard J., and Neo D. Martinez. 2000. “Simple Rules Yield Complex Food Webs.” Nature 404 (6774): 180–83. https://doi.org/10.1038/35004572.

Williams, Richard J., and Neo D. Martinez. 2004. “Limits to Trophic Levels and Omnivory in Complex Food Webs: Theory and Data.” The American Naturalist 163 (3): 458–68. https://doi.org/10.1086/381964.

Yodzis, Peter, and Kirk O. Winemiller. 1999. “In Search of Operational Trophospecies in a Tropical Aquatic Food Web.” Oikos 87 (2): 327–40. https://doi.org/10.2307/3546748.