A dynamic mechanism for prevalence of triangles in competitive networks

A dynamic mec hanism for prev alence of triangles in comp etitiv e net w orks M.N. Mo oij 1,2,* , M. Baudena 3,4,5 , A.S. v on der Heydt 2,6 , L. Miele 3,7 , and I. Kryv en 1,2,** 1 Mathematical Institute, Utrec ht Univ ersity , Utrec ht, The Netherlands 2 Cen tre for Complex Systems Studies (CCSS), Utrech t Universit y , Utrech t, The Netherlands 3 National Researc h Council of Italy , Institute of A tmospheric Sciences and Climate (CNR-ISAC), 10133 T orino, Italy 4 National Bio div ersit y F uture Cen ter (NBF C), 90133 P alermo, Italy 5 Cop ernicus Institute of Sustainable Dev elopmen t, Utrech t Universit y , Utrech t, The Netherlands 6 Institute for Marine and A tmospheric Research (IMA U), Utrech t Universit y , Utrech t, The Netherlands 7 INRAE, UR1115, Plan tes et Systèmes de culture Horticoles (PSH), Site Agroparc 84914, F rance * m.n.mo oij@uu.nl ** i.v.kryv en@uu.nl Abstract T riangles are abundan t in real-world net works but rar e in standard null mo dels for sparse graphs. Existing explanations t ypically rely on explicit triadic closure mec hanisms or geometry-based connection rules. W e prop ose an alternative h yp othesis: the frequen t app earance of triangles ma y arise naturally from the requiremen t of dynamic stability that main tains co existence of sp ecies in Lotk a–V olterra systems with comp etitive interactions. T o ev aluate this idea, we show that, across all p ossible in teraction graphs, co existence is guaran teed whenev er the coupling strength is b elow the recipro cal of the graph’s maximum degree. W e also show that co existence can p ersist up to a critical coupling strength of 1, which lea ves a large gap that is unexplained b y the graph degrees alone. These low er and upp er b ounds are ac hieved for star and complete graphs resp ectively . T o inv estigate what structural prop erties of the interaction graph control the critical coupling within the gap, w e optimise netw orks algorithmically while keeping the degree sequence fixed . W e find that net works supp orting stronger interaction strengths consis- ten tly exhibit higher clustering coefficients in sev eral netw ork mo dels. Moreo ver, in real-w orld grassland plan t netw orks, w e observ e higher clustering and stronger stability than those ex- p ected from a configuration model with the same degree sequence. Our result suggests that triangles, and clustering in general, may emerge as a structural signature of stabilising comp e- tition. Keyw ords: T riadic closure, Stabilit y in Complex Systems, Comp etitiv e Dynamics, Net work Optimisation. 1 1 In tro duction In sufficiently sparse net work mo dels where edges are placed indep enden tly at random, triangles b ecome increasingly rare as the net work size grows [1]. This contrasts with observ ations in man y large real-world net works, which hav e more triangles than their corresponding n ull models with the same degree sequence, for example, see the discussion in [2, 3] and the references therein. While triangles are una v oidable in dense net works due to the threshold and saturation effects [1], in sparse net works their abundance is typically attributed to t wo main mec hanisms. In so cial netw orks, this phenomenon is commonly explained b y triadic closur e – the increased probabilit y that t wo v ertices become connected when they ha v e a common neigh b our. In particular, this is reflected in mo dels with strong lo cal clustering, suc h as the W atts–Strogatz small-w orld mo del [2]. In geometric net works, where v ertices ha v e coordinates in a laten t metric space, lo cal clustering arises naturally from the triangle ine quality : if tw o vertices are b oth close in latent space to a third one, they are also lik ely to b e close to each other and hence to induce a (top ological) triangle [3]. Finally , it is worth mentioning that simply conditioning a generic random graph mo del, suc h as the Erdős–Rén yi random graph or exp onential random graph mo dels, to contain many triangles typically pro duces highly regular and hence unrealistic structures as a consequence of measure degeneracy [4, 5]. Man y net works that feature clustering, how ev er, do not fall under either of the ab ov e-men tioned categories. In this pap er w e give a p ossible explanation for the surplus of triangles using asymp- totic stabilit y theory for dynamical systems with comp etitive interactions. Namely , we consider a second-order system of ODEs gov erned by the Lotk a-V olterra (L V) equations on a netw ork with symmetric negative in teractions. The in teraction strength, referred to as the c omp etitive pr essur e , and the carrying capacity are identical for all species. T o quan tify stabilit y , we deter- mine the maxim um comp etitive pressure (in absolute v alue) that the system can sustain while still admitting a coexistence equilibrium, that is, an equilibrium in which all species main tain strictly positive abundances in the large-time limit. W e iden tify the maxim um comp etitiv e pressure using a bifurcation analysis and call it the criti c al c oupling of the system, with larger v alues indicating a more robust co existence equilibrium. The relation betw een the critical coupling and the structure of the netw ork has b een a topic of extensive researc h. Rob ert Ma y argued that the critical coupling of systems with a w eighted Erdős–Rén yi in teraction matrix tends to zero as the system size grows [6]. Gao et al. observ ed that a bifurcation at the critical coupling in large dimensional systems with random in teractions is captured b y a one-dimensional ‘effective’ manifold [7]. In our previous w orks, w e connected the critical coupling to Katz cen trality [8] and sho wed that this quantit y is b ounded a wa y from zero when the maxim um degree (or the maximum total weigh t around a v ertex) remain finite [9]. In th is work, we dev elop an extremal graph theory argument to show that the critical cou- pling can tak e v alues only in the w ell-defined in terv al [∆ − 1 , 1] , where ∆ is the maxim um degree of the netw ork. W e then fo cus on several random graph mo dels and show that maximising the clustering while keeping the degrees fixed also tends to increase the critical coupling in net- w orks sampled from these mo dels, and the reverse also holds: optimising for critical coupling enhances clustering. W e then turn to a real-w orld netw ork inferred from herbaceous plant comm unities across Northern Eurasia [10] and show that b oth its clustering and its critical coupling are larger than what would b e exp ected in the maxim um-entrop y n ull mo del with an iden tical degree sequence. T ak en together, these results suggest that the abundance of triangles in net works exceeding n ull-mo del predictions ma y reflect a selectiv e bias to ward systems that are more robust in the sense of asymptotic stabilit y . 2 0.00 0.25 0.50 0.75 1.00 τ 0.00 0.25 0.50 0.75 1.00 x min complete path star Figure 1: Least abundant sp ecies vs. coupling strength. The dep endence of the minim um-abundance sp ecies v alue x min at equilibrium on the coupling strength τ is sho wn for three graphs with 6 no des: complete, path, and star. Abrupt jumps in the equilibrium branc h app ear as discontin uities, o ccurring at the critical coupling where the equilibrium loses feasibilit y . 2 Results W e consider a comp etitive L V system on a fixed, but p ossibly random, undirected net work with a symmetric adjacency matrix A , d x i d t = x i (1 − x i ) − τ n X j =1 A ij x i x j , (1) where x ∈ R n is the vector of sp ecies abundances and the in teraction strength τ ≥ 0 represen ts the comp etitive pressure. In the absence of interactions b etw een the no des, meaning τ = 0 , the system reduces to n indep enden t logistic equations, each with a carrying capacity of 1. In this case, the system will con v erge globally to an equilibrium state x ∗ , where x ∗ i = 1 for all i = 1 , 2 , . . . , n . The in tro duction of no de in teractions significan tly alters the dynamics, with the system b ehaviour no w dep ending on the underlying netw ork structure. The stabilit y of system (1) is gov erned by its Jacobian, J ( x ) ik = 1 − 2 x i − τ n X j =1 A ij x j ! 1 { i = k } − τ x i 1 { A ik =1 } . (2) Here, 1 {·} denotes the indicator function (equal to 1 if the condition holds and 0 otherwise); thus 1 { i = k } selects the diagonal terms, while 1 { A ik =1 } selects the off-diagonal entries corresp onding to edges in the netw ork. If it exists, the equilibrium fixed p oin t with x ∗ i > 0 for all i = 1 , 2 , . . . , n is unique, globally stable [11], and given by x ∗ ( τ ) = ( I + τ A ) − 1 1 . (3) W e refer to x ∗ ( τ ) as a branch of feasible points. Along this branch, the Jacobian simplifies to J ( x ∗ ) ik = − x i  1 { i = k } + τ 1 { A ik =1 }  . (4) 3 0 ∆ − 1 τ c 1 guaran teed co existence guaran teed extinction G − G + critical coupling smaller clustering larger clustering Figure 2: Bounds for τ c , the maximum coupling p ermitting co existence. A cross all graphs G , the lo wer and upp er b ounds are attained at G − and G + , respectively , where the definitions of G − and G + dep end on the constraint imposed on G – either a fixed num b er of no des or a fixed maximum degree, as explained in Theorem 1. Consider τ contin uously increasing on the in terv al [0 , τ c ) . This giv es rise to tw o distinct types of bifurcation p oints that dep end on the net w ork structure [9]: (1) A b oundary transcritical bifurcation occurs when the branc h of stable feasible points x ∗ ( τ ) collides with one (or several) of the inv arian t hyperplanes H i := { x ∈ R n : x i = 0 } . (2) A saddle node or pitchfork bifurcation o ccurs when the feasible equilibrium branc h x ∗ ( τ ) loses stabilit y as the leading eigenv alue of the Jacobian crosses zero. By definition, this ev ent happ ens when the branch is in the interior of the domain. In b oth scenarios the feasible fixed p oin t loses stability , and the corresponding equilibrium branc h b ecomes confined to one of the b oundary inv ariant hyperplanes H i (or to an in tersection of several suc h hyperplanes). This even t is in terpreted as the extinction of one (or sev eral) sp ecies. W e define the critic al c oupling τ c as the smallest v alue of τ > 0 for which at least one comp onen t of the fixed p oin t x ∗ i ( τ ) is zero, which corresp onds to at least one extinction. τ c = inf τ > 0 { x ∗ i ( τ ) = 0 , i ∈ { 1 , . . . , n }} . (5) As illustrated in Figure 1, the critical couplings correspond to discontin uities at the bifurcation p oin ts. In Section 2.1 w e derive low er and upp er b ounds for the critical coupling. F urthermore, w e sho w in Section 2.2 using optimisation algorithms, that the critical coupling and mean clustering are strongly dep enden t in sev eral random graph mo dels. Finally , Section 2.3 sho ws a similar observ ation in real-w orld comp etition netw orks. 2.1 Univ ersal b ounds for critical coupling In this section, we deriv e lo w er and upp er b ounds for the critical coupling and show that τ c ∈ [∆ − 1 , 1] for any graph with maxim um degree ∆ . In particular, this shows that the lo wer b ound on τ c is independent of the net w ork size. Consequen tly , arbitrarily large biological systems can b e stable provided their maxim um degree remains b ounded. W e further show that for a fixed n umber of no des n , the least robust structure is the star graph S n , while the most robust is the complete graph K n . These results are formalized in the follo wing main theorem, whic h is illustrated in Figure 2. 4 Theorem 1. F or every c onne cte d gr aph G = ( V , E ) on n vertic es with | E | > 0 , τ c ( G ) ∈  ( n − 1) − 1 , 1  , (6) wher e the lower b ound is uniquely attaine d by the star gr aph S n and the upp er b ound by the c omplete gr aph K n . F or every finite gr aph G with maximum de gr e e ∆ > 0 , τ c ( G ) ∈  ∆ − 1 , 1  , (7) wher e the lower b ound is attaine d whenever G c ontains a c omp onent isomorphic to S ∆+1 , while the upp er b ound is attaine d by gr aphs whose c onne cte d c omp onents ar e al l c omplete gr aphs. The proof of this theorem is an immediate consequence of Lemma 1, Lemma 2, and Lemma 3, whic h are form ulated in the rest of this section. Figure 2 illustrates the statements of Theorem 1. Lemma 1. L et G b e a simple gr aph on n vertic es with maximum de gr e e ∆ , and let τ c ( G ) denote the critic al c oupling. Then ∆ − 1 ≤ τ c ( G ) ≤ − λ min ( A ) − 1 . Pr o of. Let A b e the adjacency matrix of G . The feasible branch is given by x ∗ ( τ ) = ( I + τ A ) − 1 1 , whic h is w ell-defined for τ < ∆ − 1 ≤ ρ ( A ) − 1 . F or τ = 0 this branc h is stable. As τ increases, either a saddle-no de or pitc h fork bifurcation o ccurs when the leading eigenv alue of the Jacobian J ( x ∗ ) = − diag ( x ∗ )( I + τ A ) crosses zero, whic h happ ens when τ = − λ min ( A ) − 1 . Hence no such bifurcation o ccurs for τ < ∆ − 1 ≤ ρ ( A ) − 1 ≤ − λ min ( A ) − 1 . It remains to rule out a transcritical bifurcation for τ < ∆ − 1 . W e write the fixed p oint comp onen twise, x ∗ i = 1 − τ X j A ij x ∗ j ≤ 1 , and therefore x ∗ i ≥ 1 − τ deg( i ) ≥ 1 − τ ∆ . Th us x ∗ i ( τ ) ≥ 1 − τ ∆ ≥ 0 , for all i = 1 , . . . , n, and so x ∗ ( τ ) stays aw a y from the in v ariant hyperplanes for these v alues of τ . Lemma 2. Among al l c onne cte d gr aphs on n vertic es, the star gr aph S n uniquely minimises the critic al c oupling, with τ c ( S n ) = ( n − 1) − 1 . (8) Pr o of. F or the star graph, the smallest comp onent of x ∗ o ccurs at the central no de. Expanding x ∗ ( τ ) into its series, x ∗ k ( τ ) = ∞ X i =0 ( − 1) i τ i W i ( k ) , 5 where W i ( k ) denotes the num b er of w alks of length i starting at no de k . F or the cen tral v ertex, denoted 1 , the equilibrium v alue is x ∗ 1 ( τ ) = ∞ X i =0 ( − 1) i τ i ( n − 1) ⌊ ( i +1) / 2 ⌋ . Separating ev en and o dd terms giv es geometric sums x ∗ 1 ( τ ) = ∞ X j =0 τ 2 j ( n − 1) j − ∞ X j =0 τ 2 j +1 ( n − 1) j +1 = 1 1 − ( n − 1) τ 2 − ( n − 1) τ 1 − ( n − 1) τ 2 = 1 − ( n − 1) τ 1 − ( n − 1) τ 2 . The smallest p ositive ro ot is therefore τ c ( S n ) = ( n − 1) − 1 . T o prov e uniqueness, supp ose that for some connected graph G w e hav e x ∗ i = 0 at τ = ( n − 1) − 1 . Then n X j =1 A ij x ∗ j = τ − 1 = n − 1 . Since eac h x ∗ j ≤ 1 , it follo ws that x ∗ j = 1 and deg ( i ) = n − 1 for all neigh b ours j of i . Th us S n cen tred at i is a subgraph of G . F or any neighbour j of i , 1 + τ n X k =1 A j k x ∗ k = 1 , so n X k =1 A j k x ∗ k = 0 . Since x ∗ k ≥ 0 b y contin uit y of x ∗ ( τ ) , every neighbour of j must hav e v alue 0 . Hence j is only connected to i , and therefore G = S n . Corollary 1. F or any simple gr aph G with maximum de gr e e ∆ , τ c ( G ) ≥ ∆ − 1 . If one of the c onne cte d c omp onents of G is a star S ∆+1 , then e quality holds. Pr o of. The lo wer b ound follo ws from Lemma 1. Recall that for a disconnected graph, the critical coupling equals the minimum of the critical couplings of its connected comp onen ts. If G has connected comp onent C ∼ = S ∆+1 , then τ c ( G ) ≤ τ c ( S ∆+1 ) = ∆ − 1 . Com bined with the general lo wer b ound τ c ( G ) ≥ ∆ − 1 , w e conclude that τ c ( G ) = ∆ − 1 . Remark 1. Note that the sufficient c ondition in the ab ove c or ol lary is not ne c essary. Inde e d, c onsider the double star on 7 vertic es: two hubs i, ℓ of de gr e e 3 sharing a c ommon neighb our j 1 , with two private le aves e ach (se e Figur e 3). Then ∆ = 3 and x ∗ i ( τ ) = x ∗ ℓ ( τ ) = 1 − 3 τ 1 − 4 τ 2 , which vanishes at τ = 1 3 = ∆ − 1 . Thus τ c = ∆ − 1 , even though no c omp onent in the gr aph is a star S 4 . 6 i ℓ j 1 j 2 j 3 j 4 j 5 Figure 3: The double star with ∆ = 3 . Shaded nodes are the h ubs where x ∗ v anishes first. Lemma 3. L et G b e a simple gr aph on n vertic es, and let τ c ( G ) denote the critic al c oupling. Then τ c ( G ) ≤ 1 . (9) If G is c onne cte d, the e quality is attaine d if and only if G ∼ = K n . If G is not c onne cte d, the e quality is attaine d if and only if every c onne cte d c omp onent of G is a c omplete gr aph K n i with n i ≥ 2 . Pr o of. Let G b e connected and d -regular, then x ∗ i = (1 + τ d ) − 1 > 0 , for all i = 1 , . . . , n, for all τ > 0 , so no transcritical bifurcation o ccurs. In particular, let G b e a complete graph K n , then λ min ( A ) = − 1 , hence τ c ( K n ) = − 1 λ min ( A ) = 1 . Ho wev er, K n uniquely maximises λ min ( A ) among all connected simple graphs with n v er- tices [12]. Since τ c ( G ) ≤ − λ min ( A ) − 1 , as sho wn in Lemma 1, it follows that τ c ( G ) ≤ τ c ( K n ) = 1 . If G is not connected, the same reasoning applies comp onent-wise. 2.2 Net w ork optimisation T o quantify triangle density , we use the clustering co efficien ts C ( G ) = 1 n n X i =1 C i , C i = 2 T i d i ( d i − 1) , d i ≥ 2 , (10) where T i is the n umber of triangles inciden t to no de i , so that C i measures the fraction of p ossible triangles around no de i that are realized (and C i = 0 when d i < 2 ). The degree sequence of G is denoted by d G = ( d 1 , . . . , d n ) . Theorem 1 suggests that complete graphs, for whic h C ( G ) = 1 , are also the most robust in terms of ha ving maximal τ c . T o in v estigate potential dep endencies betw een the critical coupling τ c and the mean clustering C , w e optimise the net work structure with respect to these quan tities. Sp ecifically , among all graphs with degree sequence d G , w e identify those that attain the smallest and largest p ossible mean clustering by solving arg min G ′ C ( G ′ ) and arg max G ′ C ( G ′ ) s.t. d G ′ = d G . (11) W e will find these optimal netw ork structures b y setting up the double-edge-sw ap Mark o v c hain that rewires edges ( i, j ) and ( k , l ) in t wo distinct w ays: 1) p erp endicular rewiring, to ( i, k ) and ( j, l ) , and 2) crosswise rewiring to ( i, l ) and ( j, k ) . These rewiring operations preserv e the degree sequence of the net work, ensuring that eac h no de retains the same num b er of neigh b ours (Figure 4). More imp ortantly , the double-edge-swap Marko v c hain is kno wn to b e irreducible 7 crossed switched 1 2 3 4 1 2 3 4 1 2 3 4 Figure 4: Illustration of the degree-preserving edge-sw ap (rewiring) op eration. T w o edges connecting no des (1 , 2) and (3 , 4) are replaced by either crossed or switched configurations, conserving the degree of each no de while altering net w ork top ology . [13], which guaran tees that an y netw ork with the same degree sequence is reac hable. In practice, the optimisation process may also sp end a long time in v arious metastable states. Figure 5 illustrates the optimisation pro cedure for a sp ecific netw ork that can b e divided into disjoin t cliques. W e observ e that optimisation of the mean clustering of the net w ork increases the critical coupling, which aligns with the statement of Lemma 3. -15 -10 -5 0 5 10 15 Rewiring number 0.20 0.40 0.60 0.80 1.00 C r i t i c a l c o u p l i n g , c b c Mean clustering 0.00 0.20 0.40 0.60 0.80 1.00 Mean clustering a Decreased Increased Rewiring number -15 -10 -5 0 5 10 15 Figure 5: Degree-preserving rewiring that increases clustering also increases the critical coupling. Single net work with n = 14 no des and m = 20 edges. The rewiring n umber coun ts successive degree-preserving edge sw aps relativ e to the initial netw ork (0); p ositiv e v alues increase the mean clustering co efficien t and negativ e v alues decrease it. a , Representativ e net works along the rewiring tra jectory . The dashed v ertical line marks the initial net work. Purple (left) indicates reduced clustering and green (righ t) increased clustering. b , Critical coupling τ c (left axis) and mean clustering co efficient (right axis) versus rewiring num b er. T o examine the full range of clustering and critical coupling, w e apply the rewiring pro- cedure to several random netw ork mo dels. Sp ecifically , we consider random regular netw orks (uniform degree), configuration-mo del net works with degree sequences sampled from a Poisson 8 0.0 0.2 0.30 0.32 0.34 a Regular 0.0 0.2 0.18 0.19 0.19 b Poisson 0.0 0.2 0.24 0.25 0.26 c Geometric 0.00 0.25 0.25 0.26 d W atts Strogatz 0.2 0.4 0.03 0.04 e Barabási Albert 0.0 0.1 0.28 0.30 0.33 0.35 f 0.0 0.1 0.20 0.25 0.30 0.35 g 0.0 0.2 0.25 0.30 0.35 h 0.00 0.25 0.25 0.30 0.35 i 0.1 0.2 0.10 0.20 j C r i t i c a l c o u p l i n g , c Clustering coef ficient Clustering optimization Coupling optimization Data (clustering opt., ±3 SE) Data (coupling opt., ±3 SE) Linear fit Baseline (original network) Figure 6: The critical coupling shows p ositive correlation with th e mean clustering co efficien t. Each net work is rewired to maximise or minimise the critical coupling τ c or the clustering C (1000 rewirings b oth for optimisation and minimisation). F or each netw ork model w e rep ort the mean clustering co efficient and critical coupling across 100 netw orks. All net works ha ve size n = 100 and mean degree ⟨ k ⟩ ≈ 4 . F or every netw ork t yp e, all configurations share the same degree sequence. a,f , regular ( d = 4 ); b,g , Erdős–Rényi ( p = 4 /n ); c,h , geometric ( r = 2( π ( n − 1)) − 1 / 2 ); d,i , W atts–Strogatz ( d = 4 , p = 0 . 1 ); e,j , Barabási–Alb ert ( m = 2 ). Red lines show linear fits; error bars are ± 3 s.e.m.; purple stars indicate the initial, unoptimised net works. distribution, geometric random net works, and net works generated b y the W atts–Strogatz and Barabási–Alb ert mo dels. All mo dels are constructed to ha ve the same mean degree. F or eac h mo del, w e generate 100 net w orks with size 100. Within eac h mo del, all 100 realizations share the same degree sequence but differ in structure. F or ev ery net work, w e either optimise the mean clustering or the critical coupling of the net w ork and record the measured v alues. W e observ e that the mean clustering coefficient and critical coupling are p ositively correlated across all net work models (Fig. 6). This relationship is strongest in random regular netw orks and w eak est in Barabási–Albert net works, where clustering influences critical coupling only at lo w v alues. V ariance p eaks in Erdős–Rén yi netw orks and remains small in random regular net works. Note that the optimisation pro cedure is imp ortant, as sho wn by the difference b et ween obtained v alues in the first and second ro w. 2.3 Real-w orld net works W e apply our framew ork to real-w orld net works represen ting natural grassland ecosystems, where interspecific comp etition is primarily driven b y resource limitation. W e infer these net works from the dataset compiled b y Sc heifes et al. [10], which com bines v egetation-plot comp osition with measurements of n utrien t concen trations — nitrogen (N) and phosphorus (P) — in the plants, in herbaceous plant comm unities across Northern Eurasia. W e consider 872 plots classified into habitat types according to the EUNIS classification scheme [14], sum- marized in T able 1. F or eac h habitat, the set of all sp ecies observ ed across its plots defines the corresp onding "sp ecies p o ol", from whic h we sample to construct competitive net works. Sp ecifically , we dra w random subsets of 100 species from eac h species p o ol and use them as the node set of the net work. P airwise comp etition co efficients are computed from the niche 9 o verlap of species-sp ecific resource-use strategies [15]. W e retain the strongest 10% of edges. More details are pro vided in the Metho ds section. T o test whether the resulting netw ork struc- ture can b e explained by degree heterogeneity alone, we compare each realized netw ork to a configuration-mo del ensemble that preserv es its degree sequence. F or eac h environmen t, we generate 500 configuration-mo del realizations as simple net works (i.e., no self-loops and no m ulti-edges) with the same degree sequence as the realized net work. W e then compare the realized clustering co efficien t difference ∆ C and the critical coupling difference ∆ τ c (Figure 7a) against the corresp onding distributions under the null ensemble. In all environmen ts, ∆ C > 0 and ∆ τ c > 0 for all M = 500 configuration-mo del realizations. T able 1: Habitat classification acron ym, group and t yp e, with the n umber of sp ecies recorded p er each habitat. V ar. Habitat group Habitat t yp e # species R1 Grasslands and lands dominated b y forbs, mosses or lic hens Dry grasslands 128 R2 Mesic grasslands 230 R3 Seasonally w et and wet grasslands 303 N1 Coastal habitats Coastal dunes and sandy shores 158 Q2 W etlands V alley mires, p o or fens and transition mires 155 Q5 Heloph yte b eds 192 S9 Heathlands, scrub and tundra Riv erine and fen scrub 151 0.42 0.44 0.46 0.48 0.50 0.52 C 0.040 0.045 0.050 0.055 0.060 0.065 c a N1 Q2 Q5 R1 R2 R3 S9 b 0.0 0.5 1.0 Interaction strength 0 400 800 1200 1600 Frequency c Figure 7: Construction and prop erties of realized interaction net w orks. a, Differences b et ween real net works and degree-preserving configuration-model con trols. Horizon tal axis ∆ C = C real − C conf and v ertical axis ∆ τ c = τ c, real − τ c, conf . C real and τ c, real are measured on the real net work. C conf and τ c, conf are the corresp onding configuration-mo del quantities av eraged o ver 500 realizations. Symbols denote en vironments as defined in T able 1. b, Example realized net work of size n = 100 after retaining the strongest 10% of edges, where no de size and colour are proportional to degree. c, Distribution of pairwise in teraction strengths. The opaque region indicates in teractions set to zero corresp onding to absen t edges in the realized net work; 90% of edges are remo v ed. 10 3 Discussion T riadic closure, quan tified by the clustering coefficient C , stabilises comp etitiv e netw orks. W e measure stabilit y by the critical coupling, i.e. the comp etition strength at whic h exclusions first app ear. T o separate top ology from degree effects, we compare against random-netw ork con trols that preserv e the degree sequence while changing which pairs are link ed. Under this constraint the direction is consisten t across all simulations we report: higher clustering corresp onds to larger critical coupling. Real-world grassland net works follow the same trend, b eing more clustered and showing higher critical coupling than their degree-matched controls. This suggests that adaptation migh t b e o ccurring in these netw orks. Our analysis has b een done for sparse and comp etitive netw orks, whic h are perv asive in nature, e.g. [16]. Extending the analysis to non-sparse netw orks, where sp ecies interaction is ubiquitous, could b e an interesting follow-up. Moreo v er, it still remains an op en question if other types of systems, i.e. m utualistic, facilitativ e, or trophic, also sho w the same b ehaviour. Analysing additional real-w orld net works, for example in gut microbiome [17, 18], can further strengthen our w ork. In our analysis, w e allo w for the net works to b ecome disconnected. This fragmen tation can raise the critical coupling by con taining extinction ev ents in sub-netw orks, rather than allo wing cascades to propagate through the full net w ork. That mec hanism is ecologically plausible: fragmen ted habitats and limited disp ersal w ould lo calise comp etitive effects in m uc h the same w ay . Finally , let us mention that our results are p otentially relev an t to v arious disciplines outside of ecology . In economics, the Lotk a-V olterra system is often considered. Examples are mo delling economies and banking systems [19, 20]. Also in computer science this system is used, for example in congestion control [21, 22] or p erforming com binatorial optimisation [8]. Ov erall, w e conclude that stability is not only dictated by no de degrees. The arrangement of competitive links matters, and triadic closure in particular shifts the threshold at whic h comp etitiv e exclusion sets in. 4 Metho ds Optimisation routine. W e optimise by degree-preserving edge rewiring with a Metrop olis- Hastings algorithm. Let O ( G ) denote the ob jective (mean clustering or critical coupling when increasing it, or its negative when decreasing it), and let δ = O ( G ′ ) − O ( G ) b e the change from a prop osed rewiring G → G ′ . W e accept prop osals with probability P acc = ( 1 , δ ≥ 0 , exp( δ /T ) , δ < 0 , (12) where T > 0 is a temp erature controlling acceptance of do wnhill mov es (larger T implies more exploration). This sto chastic acceptance helps escap e lo cal minima. F or eac h netw ork and for each direction (increase vs. decrease of mean clustering) w e p er- form man y rewirings. At each step we sample candidate edge pairs uniformly at random and reject an y mo ve that would create self-lo ops or m ulti-edges, thereby preserving the degree se- quence. F or every netw ork t yp e w e run multiple independent optimisations that share the same degree sequence, obtained by generating one instance of the mo del and re-randomising edges while preserving that sequence. 11 Calculation of the critical coupling. T o determine the critical coupling, we follo w the approac h of [8], and consider the iterative up date τ k +1 = τ k − n Q i =1 x i n P j =1 ( − M − 1 AM − 1 ) Q i  = j x i , (13) where M = I + τ A , leading to the bifurcation v alue τ c = lim k →∞ τ k . The equilibrium p oint x ∗ > 0 can be obtained b y inv erting the matrix M . This in verse exists for all coupling strengths τ smaller than | λ min ( A ) | − 1 , where λ min ( A ) is the smallest eigenv alue of A . When no transcritical bifurcation o ccurs first, the in verse ceases to exist at τ = | λ min ( A ) | − 1 and the system undergo es a pitchfork bifurcation. Grassland comp etitiv e in teraction net works Our aim is to obtain a data-based approx- imation of comp etitive structure within each en vironment using the nutrien t-limitation infor- mation av ailable in the dataset. The common ecological assumption is that sp ecies are more lik ely to compete strongly when their en vironmen tal preferences ("niches") are more similar, i.e. their niches o verlap more [23, 24]. T o implement this, w e defined their nic hes based on the n utrient concen trations measured in the plots where they o ccurred. If they tend to o ccur under similar n utrien t concen tration conditions, their comp etition is stronger. In particular, we use the ratio of nitrogen and phosphorus N : P as the relev ant nic he to define comp etition, as it has b een identified as a key driv er of the dynamics of these grasslands [25, 26, 27]. F or each sp ecies i , we collect the set of N : P v alues from plots in which the sp ecies is present, where presence is defined as abundance > 0 . W e then calculate the mean and standard deviation of these v alues, denoted ( µ i , σ i ) . Finally , w e assume that the niche of the sp ecies i along the N : P axis is a normal distribution f i ( x ) with mean µ i and standard deviation σ i where x := N : P . F ollowing a common approach in ecology [15, 28], pairwise interaction strength is given b y the nic he ov erlap, quantified by the Piank a index [29]. W e retain only the strongest inferred in teractions, remo ving the weak est 90% of off-diagonal o verlap v alues and keeping the strongest 10% (Figure 7a). W e p erformed a sensitivit y analysis on the exact p ercentage of edges remov ed, see b elow. After thresholding, retained ov erlaps are treated as edges and remo v ed o verlaps as absen t edges. The realized net work is represen ted b y an unw eigh ted adjacency matrix A obtained b y binarising the thresholded o verlaps; self- in teractions are excluded ( A ii = 0 ). An example realized net work with n = 100 is sho wn in Figure 7b. Sensitivit y analysis of real-w orld netw orks. T o assess if our results w ere sensitiv e to the exact p ercen tage of remov ed edges, w e repeated the analysis for the p ercentages 75% , 80% , 85% , 90% , 95% for three habitats (R1, Q5 and S9). Each sampled net work has 100 nodes; for eac h habi- tat/p ercen tage w e generated 50 sampled netw orks and, for each, 100 configuration-mo del net- w orks with matched degree sequence. W e observe that clustering and the critical coupling are closer to the v alues predicted by the configuration mo del. How ev er, they remain distinctly higher o verall (Fig. 8). 5 Co de a v ailabilit y The co de for this pap er is a v ailable at http s:/ /git hub .co m /Ni ekM o oij /An- a dapt ati on- me cha nism - reg ula ting- triad ic- c los u re . 12 0.40 0.45 0.50 0.55 C 0.03 0.04 0.05 0.06 0.07 c Habitat R1 Q5 S9 Percentage zero 75% 80% 85% 90% 95% Figure 8: Sensitivit y of the analysis to the p ercentage of edges that are remov ed across habitats. Poin ts sho w p ercentage means for each habitat (blue/orange/green), with error bars indicating ± 1 standard error. Lines connect percentages within eac h habitat; mark er shap es denote the p ercentage of edges remov ed. Each sampled net work has 100 nodes; for each habitat/p ercen tage we generated 50 sampled netw orks and, for each, 100 configuration-mo del net works with matched degree sequence. 6 Data a v ailabilit y No data is collected for this manuscript. References [1] Bollobás, B. R andom Gr aphs (Cam bridge Universit y Press, Cambridge, 2001), 2 edn. [2] W atts, D. J. & Strogatz, S. H. Collectiv e dynamics of ‘small-world’ net works. Natur e 393 , 440–442 (1998). [3] Boguna, M. et al. Netw ork geometry . Natur e R eviews Physics 3 , 114–135 (2021). [4] Chakrab ort y , S., v an der Hofstad, R. & den Hollander, F. Sparse random graphs with man y triangles. Indagationes Mathematic ae (2026). [5] Chatterjee, S. & Diaconis, P . Estimating and understanding exponential random graph mo dels. The Annals of Statistics 41 , 2428–2461 (2013). [6] Ma y , R. M. Will a large complex system be stable? Natur e 238 , 413–414 (1972). [7] Gao, J., Barzel, B. & Barabási, A.-L. Univ ersal resilience patterns in complex netw orks. Natur e 530 , 307–312 (2016). [8] Mo oij, M. N. & Kryven, I. Finding large indep endent sets in net w orks using comp etitiv e dynamics. Communic ations Physics 8 , 235 (2025). 13 [9] Mo oij, M. N., Baudena, M., v on der Heydt, A. S. & Kryven, I. Stable co existence in indefinitely large systems of competing sp ecies. Pr o c e e dings of the R oyal So ciety A: Math- ematic al, Physic al and Engine ering Scienc es 480 , 20240290 (2024). [10] Sc heifes, D. J. P . et al. The plan t ro ot economics space in relation to nutrien t limitation in Eurasian herbaceous plant comm unities. Ec olo gy L etters 27 , e14402 (2024). [11] Goh, B. S. Global stabilit y in man y-sp ecies systems. The Americ an Natur alist 111 , 135– 143 (1977). [12] Brou wer, A. E. & Haemers, W. H. Sp e ctr a of Gr aphs (Springer, New Y ork, 2012). [13] Greenhill, C. The switc h Mark ov chain for sampling irregular graphs. Combinatorics, Pr ob ability and Computing 24 , 556–573 (2015). [14] Da vies, C. E., Moss, D. & Hill, M. O. EUNIS habitat classification revised 2004. R ep ort to: Eur op e an envir onment agency-Eur op e an topic c entr e on natur e pr ote ction and bio diversity 127–143 (2004). [15] Haraldsson, M. & Thébault, E. Emerging nic he clustering results from b oth competition and predation. Ec olo gy L etters 26 , 1200–1211 (2023). [16] Connell, J. H. On the prev alence and relativ e imp ortance of in tersp ecific competition: Evidence from field exp eriments. The A meric an Natur alist 122 , 661–696 (1983). [17] Co yte, K. Z., Sc hluter, J. & F oster, K. R. The ecology of the microbiome: Net works, comp etition, and stability . Scienc e 350 , 663–666 (2015). [18] Camac ho-Mateu, J., Lamp o, A., Sireci, M., Muñoz, M. A. & Cuesta, J. A. Sparse sp ecies in teractions repro duce abundance correlation patterns in microbial comm unities. Pr o c e e d- ings of the National A c ademy of Scienc es 121 , e2309575121 (2024). [19] Comes, C.-A. Banking system: Three level Lotk a-V olterra mo del. Pr o c e dia Ec onomics and Financ e 3 , 251–255 (2012). [20] Moran, J. & Bouchaud, J.-P . May’s in stabilit y in large economies. Phys. R ev. E 100 , 032307 (2019). [21] An toniou, P . & Pitsillides, A. Congestion control in autonomous decentralized net works based on the Lotk a-V olterra comp etition mo del. In Alippi, C., P olycarp ou, M., Pana yiotou, C. & Ellinas, G. (eds.) Artificial Neur al Networks – ICANN 2009 , 986–996 (Springer Berlin Heidelb erg, Berlin, Heidelb erg, 2009). [22] Iguc hi, T., Hasega w a, G. & Murata, M. A new congestion control mechanism of tcp with inline netw ork measurement. In Kim, C. (ed.) Information Networking. Conver genc e in Br o adb and and Mobile Networking , 109–121 (Springer Berlin Heidelb erg, Berlin, Heidel- b erg, 2005). [23] W andrag, E. M., Catford, J. A. & Duncan, R. P . Quantifying niche av ailability , nic he ov er- lap and comp etition for recruitment sites in plant p opulations without explicit knowledge of nic he axes. Journal of Ec olo gy 107 , 1791–1803 (2019). [24] Chesson, P . Mechanisms of maintenance of sp ecies div ersity . Annual r eview of Ec olo gy and Systematics 31 , 343–366 (2000). 14 [25] F ay , P . A. et al. Grassland pro ductivity limited b y multiple n utrients. Natur e Plants 1 , 1–5 (2015). [26] Güsew ell, S. N: P ratios in terrestrial plants: v ariation and functional significance. New Phytolo gist 164 , 243–266 (2004). [27] Olde V en terink, H., W assen, M. J., V erkro ost, A. & De Ruiter, P . Sp ecies ric hness– pro ductivit y patterns differ b et ween n-, p-, and k-limited wetlands. Ec olo gy 84 , 2191–2199 (2003). [28] Mason, N. W., Mouillot, D., Lee, W. G. & Wilson, J. B. F unctional richness, functional ev enness and functional div ergence: the primary comp onen ts of functional div ersity . Oikos 111 , 112–118 (2005). [29] Piank a, E. R. The structure of lizard comm unities. Annual r eview of Ec olo gy and System- atics 53–74 (1973). 7 F unding This w ork w as supp orted by the Italian National Biodiversit y F uture Cen ter (NBF C): Na- tional Recov ery and Resilience Plan (NRRP), Mission 4 Comp onent 2 Inv estment 1.4 of the Italian Ministry of Univ ersit y and Researc h; funded by the EU - NextGenerationEU (Pro ject co de CN 00000033). MB also ac kno wledges funds from the SENTINEL PRIN pro ject – Call 2022, Grant No. 2022CM4F3X. This w ork was supp orted b y the follo wing NW O programs: T alentprogramme, pro ject n umber VI.C.202.081; and researc h program VIDI, pro ject n um- b er VI.Vidi.213.108. F urthermore, the pro ject ackno wledges funds from the SENTINEL PRIN pro ject – Call 2022, Gran t No. 2022CM4F3X by th e Italian Ministry of Univ ersity and Researc h (MUR); from the DiviN-P pro ject under the 2021-2022 Bio divProtect join t call, co-funded b y the Europ ean Commission (GA No. 101052342) and MUR (CUP:B53C23001040006); and from the Italian National Bio div ersity F uture Cen ter (NBF C): National Recov ery and Resilience Plan (NRRP), Mission 4 Comp onen t 2 Inv estment 1.4 of MUR, funded by the EU - NextGenera- tionEU (CN 00000033). 8 A c kno wledgemen ts The authors are grateful to Annegreet V eken for discussion that help ed define the real-world net works and for supplying the data used to generate these net w orks. The authors are grate- ful to Jochem Hoogendijk, Mike de V ries, and Rub en Hendriks for helpful discussions. MNM gratefully ackno wledges support from Complex Systems F und, with sp ecial thanks to P eter K o eze. MB ac knowledges funds from the SENTINEL PRIN pro ject – Call 2022, Gran t No. 2022CM4F3X by the Italian Ministry of Univ ersity and Researc h (MUR); from the DiviN-P pro ject under the 2021-2022 BiodivProtect join t call, co-funded by the Eu rop ean Commission (GA No. 101052342) and MUR (CUP:B53C23001040006); and from the Italian National Bio- div ersity F uture Center (NBFC): National Reco v ery and Resilience Plan (NRRP), Mission 4 Comp onen t 2 Inv estmen t 1.4 of MUR, funded b y the EU - NextGenerationEU (CN 00000033). A vdH gratefully ackno wledges supp ort from NW O under the NWO T alentprogramme, pro ject n umber VI.C.202.081. IK gratefully ac knowledges supp ort from Netherlands Researc h Organ- isation (NW O), researc h program VIDI, pro ject n um b er VI.Vidi.213.108. 15 9 Comp eting in terests The authors declare no comp eting interests. 10 Materials & Corresp ondence Please con tact MNM for corresp ondence. 16