Clustering Compresses Attractors in Watts-Strogatz Threshold Boolean Networks

Does higher clustering shorten attractor periods? We examine whether the global clustering coef- ficient C, a direct measure of triangle density, predicts attractor lengths in synchronous, signed- threshold Boolean networks on Watts-Strogatz (WS) graphs. We generate330 directed, signed WS networks spanning sizesN = 10–100 and mean degrees ¯k = 2–10, simulate dynamics from M = 100 random initial states per graph with exact cycle detection, and summarise each graph by the average log attractor period (equivalently, the geometric mean period). Our primary anal- ysis relates this log-period summary toC while adjusting forN, ¯k, the mean directed shortest path (MSP), and including nonlinear size-degree and clustering-degree interactions. Main finding. Higher clustering robustly short- ens attractor periods: a 0.10 increase in C (C∈[0, 1]) corresponds to an≈13–15% lower ex- pected geometric mean period, and moving from C = 0.000 to C = 0.460 yields an≈50% reduc- tion, holding other properties fixed. The effect persists when the linearC term is replaced by a nonlinear function ofC, and it replicates in held-out graph instances (graphs not used to fit the model). Shorter periods are not explained by an increase in fixed points under the strict comparator ( >); rather, higher triangle density shifts mass from long cycles to medium-length cycles. Significance. In threshold-like logic, settling speed and oscillatory stability are central to com- putation and control. Our results provide a di- rect, quantitative link between triangle density and these long-run behaviours, showing that C acts as a structural lever on temporal complexity. All figures and tables are reproducible from the accompanying code, data, and analysis scripts.

Boolean networks; small-world networks; clustering coefficient; attractors; graph properties; threshold dynamics; network topology. 1 Introduction 1.1 Boolean networks and attractors Boolean networks (BNs) provide a compact lan- guage for studying how the wiring of a regulatory system shapes its long-term behaviour [1, 2, 3]. Each node holds a binary state (on/off); at dis- crete time steps all nodes update according to local logical rules that encode activation and re- pression. Despite their simplicity, such models capture essential features of biological regulation when kinetic parameters are unknown or noisy, 1 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint and they offer transparent links between struc- ture and dynamics through the lens of attractors- fixed points and periodic cycles that represent stable phenotypes or recurring activity patterns [4, 5]. Kauffman’s seminal work introduced ran- dom Boolean networks as abstract models of ge- netic regulation, proposing that stable cellular phenotypes correspond to attractors in these sim- plified networks [3]. Subsequent studies formal- ized the attractor concept and showed that net- work connectivity and logic determine whether dynamics settle into an ordered model (with many fixed points and short cycles) or a chaotic model (with few fixed points and long, complex cycles) [2, 6, 7, 8, 9, 10]. Biological networks are thought to operate near this critical balance between order and chaos, ensuring both stability and adaptabil- ity, and indeed BNs have been applied to real systems to capture such dynamics. For exam- ple, a Boolean model of the yeast cell cycle net- work reproduces its robust oscillatory progression through cell cycle states [10, 11]. 1.2 Clustering and small-world structure One structural feature that may modulate a net- work’s dynamical model is clustering. Clustering measures the prevalence of closed triangles in a network, the tendency for neighbours of a node to also be connected with one another and re- flects the density of local feedback cycles [12, 13]. Many natural networks, including biological reg- ulatory networks, exhibit significant clustering [14, 15, 16]. In small-world networks, introduced by Watts and Strogatz (1998), high clustering combines with short path lengths, enabling rapid communication while maintaining local commu- nity connection [17]. Because clustering directly counts local feedback motifs (triangles) and can be tuned in the WS models, it is a natural struc- tural axis on which to study attractor statistics [17, 18]. 1.3 Prior work and mechanistic expectations Prior work has linked small-world structure to ordered dynamics, but mostly through indirect measures. However, most studies evaluate proxies (storage, damage) rather than directly quantify- ing how the clustering coefficientC relates to attractor period length. For these reasons, we hypothesise that increasing clustering should shift a Boolean network’s dynamics toward the ordered side, yielding shorter attractor cycles on average. 1.4 Modelling framework: threshold dynam- ics on directed WS graphs We address this gap by testing the cluster- ing–period link using a biologically grounded up- date rule. Threshold logic is widely used to ap- proximate regulatory interactions because it mir- rors the sigmoidal responses and saturations com- mon in gene regulation and signalling [19, 5, 20]. Majority-like threshold rules reduce effective sen- sitivity and often promote reliable settling, a tendency also observed in networks enriched for canalising functions, where a single input can de- termine the output regardless of others [21, 9, 22]. In fact, canalizing rules, which stabilize dynam- ics and prevent chaos, are frequently included in carefully constructed biological Boolean mod- els [5, 20]. Recent theory on threshold Boolean networks further clarifies stability conditions and supports using these rules to probe structure dy- namics links [23, 24]. For these reasons, we adopt synchronous signed-threshold updates and vary network structure while keeping the rule class fixed. This design isolates the effect of clustering without confounding it with changes in logical update properties. Our setting is the directed WS model, analysed on the scale of the global clustering coefficientC (see Sec. 4). The samep can produce different tri- angledensitiesacrosssizesanddegrees, whereas C is directly comparable between instances [17, 15]. Moreover, empirical networks do not come with a rewiring knob but do have measurable cluster- ing [14, 16]. In this study we estimate how the log–mean periodY = log(¯L) varies withC while adjusting for network sizeN, average degree¯k, and mean directed shortest path (MSP); we also examine average betweenness centrality (ABC) in exploratory variants to enable fair comparisons across mixed geometries [ 15, 8]. Estimation de- 2 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint tails and diagnostics are provided in Sec. 4, and substantive findings are presented in Sec. 5. 1.5 Contributions This paper makes four contributions. (1) Di- rect measurement of clustering–period links:we quantify howC relates to the attractor periodL under synchronous signed-threshold dynamics on directed WS graphs, moving beyond proxies such as AIS and damage spreading [25, 26]. (2) Effect sizes and generalisation:we report interpretable ratios on the log period scale and assess out of sample generalisation across graph instances [27]. (3)Distributional consequences:periodshortening arises primarily from a shift of probability mass away from long cycles toward medium length cy- cles; critically, with jittered thresholds ties are almost surely absent (so> and≥yield effectively identical dynamics), and the frequency of fixed points does not increase under the strict compara- tor (>), clarifying how shorter periods emerge in practice [2, 28, 4]. (4) Reproducible detection:we provide exact attractor detection via a hash-and- backtrack scheme, with parameters and stopping rules specified for replication. The method stores each visited state in a hash table and, when a repeat is found, backtracks to the first occurrence, returning the exact cycle and its period. Com- mon alternatives include SAT/BDD-based exact enumerators and state-space sampling (approx- imate), which differ mainly in memory use and scalability [4, 29]. 2 Background This section outlines the conceptual background for how clustering shapes attractor periods in threshold Boolean networks and explains why we focus on directed, signed Watts–Strogatz (WS) graphs. 2.1 Boolean networks, attractors, and phases A classic result is that random Boolean networks exhibit distinct dynamical phases depending on their parameters. In an ordered phase, pertur- bations to the state of one node tend to die out over time; such networks have many fixed points and generally short cycles. In a chaotic phase, small perturbations amplify and propagate indef- initely, resulting in long transients and complex cycles [2]. Separating these extremes is a critical edge-of-chaos model where the network is maxi- mally sensitive yet not fully chaotic [10, 30, 31]. Kauffman [3] first noted that random networks with low connectivity often settle quickly into point attractors, whereas highly connected net- works can exhibit chaotic, lengthy oscillations. Later work formalised this via the notion of av- erage sensitivity of the Boolean update rules; if the average sensitivity is less than one (for ex- ample, due to canalisation or bias towards stable outputs), the network is ordered; if it is greater than one, the network is chaotic [ 9]. Shmule- vich and Kauffman [9] quantitatively linked sen- sitivity to Boolean logic types, and Klemm and Bornholdt [7] observed that networks can have both stable and unstable attractors depending on structure and bias. These insights explain why real genetic networks are believed to operate near an ordered model [11]. Where relevant, we use signed–threshold dynamics to ground later comparisons in a biologically common rule class. 2.2 Small-world structure, clustering, and the WS baseline Beyond node-level logic, the topology of inter- actions influences a Boolean network’s dynam- ics. Many empirical networks are neither com- pletely random nor completely regular, but in- stead have small-world structure combining high clustering with short path lengths [17, 15, 25, 26]. The Watts–Strogatz (WS) model generates small- world networks by starting from a ring lattice and randomly rewiring a fraction of edges to cre- ate shortcuts [17]. At intermediate rewiring, the network retains most local links and thus high clustering while gaining a few long-range links that shrink distances [18]. Our focus on WS complements other network families. In Erdős–Rényi graphs, clustering is low and largely a function of size; mean degree is dominant [15]. In scale-free graphs, degree hetero- geneity changes sensitivity and stability; Barabási and Albert showed that hubs shape connectivity 3 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint [32], and Aldana analysed Boolean dynamics on scale-free topology [6]. Modular organisation can enhance robustness by confining perturbations within modules [33]. To vary triangle density at nearly constant degree and path length, we use WS as the ideal model [17, 18]. In other families, clustering co-varies with other graph properties, making it hard to isolate the effect of clustering alone. Our approach is complemen- tary to motif-level studies predicting attractors from small circuit templates [34]: where motif analyses are local and combinatorial, our study is global and statistical at the network level, asking how clustering relates to typical attractor periods across many graphs. 2.3 Mechanistic lenses linking clustering to periods Although different in methods, three lines of rea- soning converge on the prediction that increasing clustering should shorten attractor periods. Information storage and redundancy.In clus- tered neighbourhoods, a node’s inputs are inter- connected, creating overlapping paths of predic- tive information; small-world structure with high clustering raises active information storage (AIS) [25]. In WS networks, nodes in highly clustered regions have higher AIS, reflecting stronger local memory. Under threshold logic, such redundancy reduces effective sensitivity, supporting shorter cycles. Damage spreading and sensitivity. Damage spreading tracks how an initial perturbation prop- agates. Ferraz and Herrmann [ 35] examined Boolean dynamics on small-world graphs and ob- served that introducing even a small fraction of shortcuts could drive a transition toward more chaotic behaviour as damage spreading increased. LuandTeuscher[ 26]demonstratedthataddinglo- cal links, thereby increasing clustering, suppresses damage spreading, whereas adding long-range links, thereby decreasing clustering, amplifies it. Lower sensitivity implies shorter transients and fewer opportunities to traverse long cycles [2]. Loopclosureand feedbackmotifs. Analyses con- nect the probability of closing onto a short cycle to local feedback and redundancy [36]. Motif- level studies in related threshold linear systems show that small feedback circuits strongly con- strain possible attractors [34]. Because clustering counts triangles that form local feedback cycles, higher clustering increases opportunities for early loop closure, shifting probability toward short cycles. Together, these perspectives offer a coherent mechanistic account consistent with network the- ory [15]: triangles increase storage and decrease sensitivity; lower sensitivity and more local feed- back increase the likelihood of short cycles; there- fore, as clustering increases, attractor periods should decrease on average. Our empirical anal- ysis in Sec. 5 tests this prediction directly by relating clustering to attractor periods under a fixed update rule. 2.4 Threshold logic, canalisation, and signed edges In regulatory modelling, Boolean update rules range from unconstrained truth tables to struc- tured, biologically informed families. Thresh- old rules are widely used in gene regulation and signalling because they capture saturating, sig- moidal responses with few parameters [19, 5, 20]. Majority-type thresholds lower average sensitivity and favour ordered dynamics with reliable settling [5]. Imposing threshold logic yields robust net- work behaviour even on random topologies [37], and robustness itself has been proposed as an evolutionary organising principle [ 10]. Canalis- ing functions further stabilise dynamics and are enriched in curated biological Boolean models; networks composed entirely of canalising rules are provably stable [5, 21]. Both threshold logic and canalisation reduce effective sensitivity, tending to shorten transients and limit attractor periods [9]. Many regulatory and signalling networks are naturally represented as directed, signed graphs, with edges encoding activation (positive) or in- hibition (negative) [29, 38]. Small-world models 4 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint admit directed variants in which undirected links are oriented or edges are rewired in a directed manner, typically preserving narrow in- and out- degree distributions while retaining short paths and local cohesion [15, 18]. Measures of clus- tering generalise from undirected transitivity to directed motifs, enabling analyses that distinguish feed-forward triads from directed 3-cycles where appropriate [15, 12]. In Boolean formulations with signed thresholds, edge signs act as positive or negative weights on inputs, offering a simple and biologically interpretable way to model acti- vation and repression; when sign information is unavailable, a common neutral baseline assumes a roughly equal mix of activating and inhibiting edges [5, 37]. These considerations motivate our choice of synchronous signed-threshold dynamics on directed WS graphs (Sec. 4). 2.5 Attractor distributions, proxies for order, and structural covariates Analyses of Boolean dynamics commonly char- acterise attractor landscapes usingdistributional summaries rather than means alone. Under a specified initial state distribution and update scheme, one considers the probability of fixed points Pr[L = 1] and the probability of short non- trivial cycles (e.g.,Pr[L≤4]), together with full distributional views such as empirical cumulative distributionfunctions(ECDFs)of L[4,29]. These summariesrevealwhetherstructuralchanges(e.g., increased clustering) primarily shift mass away from long cycles toward shorter or moderate cy- cles, complementing aggregate measures such as the mean period and, for modelling purposes, its log transform [4, 27]. Much of the literature assesses dynamical or- der viaproxies such as active information storage (AIS) or damage spreading [25, 26]. A comple- mentary line of work measures attractorperiods directly and treats the period (or its log transform; see Y in Eq. 9) as the response, relating variation inL to structural features such as clustering while holding other properties in view [15, 8]. Framing

on the period scale translates qualitative expectations about storage and sensitivity into concrete, testable statements about long-run be- haviour. Attractor lengths are influenced by multiple structural features beyond clustering, including network sizeN, average degree¯k, mean (directed) shortest path (MSP; Eq. 3), and average between- ness centrality (ABC; Eq. 5) [15, 8, 39]. Because these quantities often co-vary with clustering un- der common generative models, empirical analy- ses typically adjust for them, e.g., modellinglogL as a function ofC with N, ¯k, MSP, and ABC as covariates to isolate the partial association of clustering while comparing networks that are similar in scale, sparsity, reachability, and flow centralisation. 3 Mathematics Concepts and Definitions In this section, we introduce the fundamental con- cepts and definitions that will be used throughout this paper [12, 13, 15, 18, 3, 1, 2, 4]. 3.1 Graphs, adjacency, and degree A graphG = (V,E ) has a nonemptyvertex set V and an edge setE. In an undirected graph, E⊆{{u,v}: u,v ∈V, u̸= v}; in a directed graph (adigraph),E⊆V×V with ordered pairs (u,v ). The order ofG is N =|V|and thesize of G ism =|E|. Throughout this paper, graphs are simple unless stated otherwise (no self-loops or parallel edges) [12, 13]. An adjacency matrixis A = [aij]. For an undi- rected simple graph,aij =aji = 1 if{vi,vj}∈E and 0 otherwise. For a digraph, aij = 1 if (vi,vj)∈E and symmetry is not required. For dynamical models we use asigned adjacencywith aij∈{−1, 0, +1}to encode activation(+1) and inhibition (−1). When measuringundirectedfeatures (e.g., clus- tering), we pass to theunderlying undirected sim- ple graphGu with adjacencyU = [uij] defined by uij =    1, if aij̸= 0or aji̸= 0, 0, otherwise, uii = 0for alli, and then treatGu as a simple graph (so directions and signs are both ignored and parallel edges 5 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint are not allowed). Unless otherwise noted, any undirected network statistic is computed onGu constructed fromA as above [13, 15]. The degree of vertexvi in an undirected graph is deg(vi) = ∑ juij and equals the size of the neighborhoodN(vi) ={vj :{vi,vj}∈E}. In a digraph, thein-degreeand out-degreeare degin(vi) = ∑ j 1[aji̸= 0], degout(vi) = ∑ j 1[aij̸= 0], (1) where 1[·] denotes the indicator, equal to1 when the condition holds and0 otherwise. We write the total degree of vertex vi as degtot(vi) = degin(vi) + degout(vi) (in an undirected graph, degtot(vi) = deg(vi)), and define theaverage de- gree as ¯k = 1 N ∑N i=1 degtot(vi) (so ¯k = 2m/N for simple graphs) [12]. 3.2 Paths, distance, clustering, and path- based measures A walk of length p is a sequence(u0,...,up) with adjacent consecutive vertices; apath is a walk with no repeated vertices [12]. In a digraph, we writeu⇒v to denote that there exists a directed path fromu to v. A (simple) cycle is a closed walk of length at least3 whose vertices are all distinct except the first and last. A3-cycle (a triangle) is a cycle of length3 (in a digraph, edges of the cycle must respect orientation). The distanced(u,v ) is the length of a shortest path fromu to v (respecting edge directions in digraphs). The diameter ofG is the maximum finite distance [15]. For clustering we counttriangles (3-cycles) and connected triples(wedges) on theunderlying undi- rected simple graphGu. Aconnected triple (wedge) is a path of length2,u−v−w, inGu, withu̸=w. WritingT for the number of triangles andΛ for the number of connected triples, the global clus- tering coefficient (transitivity)C is defined as: C =    3T Λ , Λ > 0, 0, Λ = 0, with C∈[0, 1], (2) which is the fraction of wedges that are closed into triangles [15, 18, 14]. We use the meandirectedshortest path (MSP), defined as: MSP =     1 |R| ∑ (u,v)∈R d(u,v ), |R|> 0, 0, |R|= 0, (3) where R ={(u,v)∈V×V : u̸=v, u⇒v in G},(4) and d(u,v ) is the length of a shortestdirected path fromu tov. Thus,G need not be (strongly) connected; the average is taken only over ordered pairs that are reachable by a directed path (i.e., pairs inR) [15]. The average betweenness centrality is ABC = 1 N ∑ v∈V ∑ s,t∈V\{v} s̸=t σst>0 σst(v) σst , (5) where σst(v) is the number of shortest directed s→tpaths that pass throughv, andσst counts shortest directeds→tpaths froms to t [40]. 3.3 Boolean networks and attractors Given a graphG = (V,E ) with|V|= N, we associate a Boolean state to each vertex at each discrete timet∈Z≥0. The state vectoris S(t) = (s1(t),...,sN(t))∈{0, 1}N, with a chosen initial conditionS(0)∈{0, 1}N. A Boolean network onG is the deterministic map F :{0,1}N→{0,1}N, S(t)↦→F (S(t)) = ( f1(S(t)),...,fN(S(t)) ) , (6) where each coordinate functionfj :{0, 1}N→ {0, 1}depends only on the states of the in- neighbors of the vertexvj in G [1, 2]. The syn- chronous evolution is then S(t+1) =F ( S(t) ) . 6 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint An update schemespecifies when node states areupdated. Inthe synchronousscheme, allnodes update simultaneously at each integer time step according to Eq.(6). Asynchronous variants (e.g., random-order or general-asynchronous) update subsets or single nodes per step; in this study we focus on synchronous updates [2, 29]. An update rulespecifies the coordinate func- tions fj used in Eq.(6) and, together with the update scheme, determines the next stateS(t+1) from S(t) [2]. A widely used rule in signed regu- latory models is thesigned–threshold rule, sj(t+1) = 1   N∑ i=1 aijsi(t) ◦θj  , (7) where aij∈{−1, 0, +1}, θj∈R is a node-specific threshold, and the comparison symbol◦∈{>,≥ }; in this study we use the strict comparator “>” throughout (see Sec. 4) [2, 5, 37]. An attractor is a minimal set A = {a0,...,aL−1} ⊆ {0, 1}N such that F (ar) = a(r+1) modL for all r∈{0,...,L−1}. Equiva- lently,FL(ar) = ar with no smaller positive pe- riod. The periodis L∈{1, 2,...}(a fixed point is the caseL = 1). The state transition graph places a directed edgeS→F (S) on the2N states (out-degree1 per state); attractors are precisely the directed cycles (with fixed points as1-cycles) [2, 4]. From a given initial state S(0), the trajectory S(0),S (1),... eventually repeats because{0, 1}N is finite. The transient lengthfrom S(0) is T = min{t≥0 : S(t)∈A}, i.e., the first time the trajectory enters the repeat- ing cycle. When summarizing outcomes over many initial states, we report theattractor periodand its log transform. FromM simulations with periodsLr, ¯L = 1 M M∑ r=1 Lr, (8) and Y = 1 M M∑ r=1 logLr, (9) so that the geometric mean attractor period is Lgeo = exp(Y ) and a difference∆Y corresponds to a multiplicative factore∆Y on Lgeo [4, 27]. 3.4 Graph models: Watts-Strogatz and Erdős- Rényi Watts-Strogatz (WS) model. Fix an integer N≥3 and an integerκ∈{1,...,⌊(N−1)/2⌋}, and let p∈[0, 1]. Let V ={0, 1,...,N −1} and foru,v ∈V with u̸= v define the circular distance δ(u,v ) = min { |u−v|, N−|u−v| } . The ring latticeL(N,κ) has edge set E0 = { {u,v}⊂V : 1 ≤δ(u,v )≤κ } , so every vertex has undirected degreedeg(u) = 2κ (i.e., κis the ring-lattice parameter). Construct the undirected WS graphH = WS(N,κ,p) by independently rewiring each{u,v}∈E0 once: with probabilityp replace{u,v}by{u,w}, where w is drawn uniformly fromV\({u}∪NH(u)) at the time of rewiring (avoids self-loops and multi- edges); otherwise keep{u,v}. This preservesm = |E(H)|=Nκand hence the mean (undirected) degree ¯k = 2m/N = 2κ[17, 15, 18]. The directed WS graphGWS =−→WS(N,κ,p) is obtained by orienting each undirected edge ofH independently and uniformly at random: P((u,v )∈E(GWS)|{u,v}∈E(H)) = 1 2, P((v,u)∈E(GWS)|{u,v}∈E(H)) = 1 2. (10) Conditioned onH, degGWS in (v)∼Binomial ( degH(v), 1 2 ) , degGWS out (v)∼Binomial ( degH(v), 1 2 ) , (11) so E [ degGWS in (v) ⏐⏐ ⏐H ] = E [ degGWS out (v) ⏐ ⏐ ⏐H ] = 1 2 degH(v). (12) 7 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint Erdős-Rényi (ER) model.Given N ≥1 and q∈[0, 1], the undirected modelH = ER(N,q ) is defined by independent edges with P ( {u,v}∈E(H) ) =q, 0≤u<v ≤N−1, (13) so E[m] =q (N 2 ) and E[¯k] =q(N−1) [15, 18]. The directed E-R graph ER→=−→ER(N,q ) is obtained by the same orientation rule as in Eq. (10). Matching expected mean degree.For compar- isons between WS(N,κ,p) and ER(N,q ), we match expected mean degree by choosing q = 2κ N−1, (14) valid whenever2κ<N, so both undirected mod- els have the same expected mean degree¯k = 2κ (and thus the same expected total degree after orientation) [15, 18]. 4 Methods 4.1 Graph model: small-world networks and clustering We use the Watts-Strogatz (WS) small-world model to vary triangle density while keeping the network sizeN and average degree¯k controlled (Sec. 3.4, Def. 3.4; [17, 15, 18]). For comparisons with the baseline we match the expected mean degree using Eq. (14). To obtain a directed, signed network consis- tent with regulatory settings, we first generate an undirected WS graph and then orient each edge once using the random orientation rule in Eq. (10), creating exactly one directed arc per undirected edge and no reciprocalu↔v pairs. We then independently assign a sign to each arc, +1 for activation and−1 for inhibition with equal probability, yielding a signed adjacency matrix A = [aij]∈{−1, 0, +1}N×N[5, 2, 15]. We measure clustering on the underlying undi- rected simple graphGu using the global clustering coefficientC (Sec. 3, Eq. (2)). 4.2 Network generation and sampling Wevariednetworksizeover N∈{10, 20,..., 100}. For each N, we considered degree inputs d ∈ {1,...,⌊N/10⌋}. Within the WS model (Def. 3.4), we set k = max{⌊d/2⌋,1}, which yields an undirected average degree¯k = 2k ∈ {2, 4, 6, 8, 10}. Each (N,d ) specification was crossed with rewiring probabilities p ∈ {0.01, 0.05, 0.10, 0.20, 0.40, 0.60}. The total num- ber of (N,d) specifications is ∑ N∈{10,20,...,100} ⌊N 10 ⌋ = 1+2+···+10 =10·11 2 = 55, so this design yields55×6 = 330 WS graphs that span topologies from highly clustered small-world networks to nearly random graphs [18]. For baseline comparisons we also generate degree-matched Erdős-Rényi graphs (Def. 3.4) with edge probabilityq chosen via Eq.(14). The main clustering results, however, are reported for the WS ensemble. Across this parameter sweep overN, ¯k, andp (with degree-matched ER base- lines), the observed clustering onGu ranges from C≈0 (nearly random) toC≈0.63 (highly clus- tered). 4.3 Graph properties and controls For each directed, signed graph, we compute five structural descriptors (Sec. 3): the sizeN; the average degree¯k; the mean directed shortest path (MSP; Eq.(3)); the average betweenness central- ity (ABC; Eq.(5)); and the clustering coefficient C (Eq. (2)) evaluated on the underlying undi- rected simple graph Gu. These quantities are defined formally in Sec. 3; here we use them to characterise each graph and as candidate covari- ates when relating clustering to dynamical out- comes [3, 2, 6, 41, 15, 18, 5, 8, 40]. In the main regression we useN, ¯k, MSP, andC as predictors; ABC is computed but only used in exploratory checks and does not enter the primary model. As an initial exploratory step, we also compute Spearman rank correlations between the aver- age log-periodY and each structural descriptor (N, ¯k, MSP, ABC,C ) to characterise simple uni- variate associations and motivate the choice of covariates in the regression models. 8 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint 4.4 Why analyse effects on theC scale? We report effects on the global clustering coeffi- cientC rather than the WS rewiring probabilityp because (i) the samep can yield different triangle densities acrossN and ¯k, (ii)C is directly compa- rable across sizes and degrees and is measurable in empirical networks, and (iii) statements on the C scale transfer beyond a specific generator [17, 15]. 4.5 Signed-threshold Boolean dynamics and thresholds We adopt synchronous updates as a standard and tractable baseline for Boolean network analyses and focus exclusively on synchronous dynamics in this study [2, 4]. We use the signed-threshold update defined in Sec. 3, Eq.(7), with a strict “>” comparator throughout (the signed adjacency A = [aij] is as defined in Sec. 3). Node thresholds follow a majority-style rule with a small tie-breaking perturbation: θj = 1 2 N∑ i=1 aij +εj, ε j∼N ( 0, 10−4 ) i.i.d., (15) which reduces average sensitivity and promotes stable behaviour in signed regulatory networks [9, 21]. The small perturbationεj breaks exact ties without materially changing the dynamics. We keep the strict comparator “>” to avoid spu- rious fixed points at equality [5, 9]. Edge signs are assigned with equal probability to provide a neutral baseline for activation-inhibition balance when system-specific information is unavailable [37, 2]. The update rule is implemented in C++ for efficient simulation and attractor detection. 4.6 Simulation and attractor detection For each graph we runM = 100 independent simulations from random initial statesS(0)∼ Bernoulli(0.5)N. This choice balances coverage of the state space and computational cost; the Bernoulli(0.5) start is an uninformative prior that matches the sign-balanced wiring [3, 2]. The fi- nite, deterministicstatespaceguaranteeseventual periodic behaviour [2]. We detect cycles exactly using a hash-based backtracking method implemented in C++. Each full state vector is stored in a hash table for constant-time lookup; when a state repeats at time t with first occurrence att1, we report the period L =t−t1 after a short backtrack to con- firm the earliest repeating state [4]. We impose a conservative cap of108 update steps per run. If a trajectory has not repeated by the cap, the run times outand is excluded from the averages for ¯L and Y; in Eq.(8) M is the number offinished runs. In our experiments, no run timed out, so all runs were included. 4.7 Dynamical outcomes If run r converges to a fixed point, its period is Lr = 1; if it enters anℓ-cycle, its period is Lr =ℓ>1. For each graph we therefore obtain M periods L1,...,LM. Our primary summary is theaverage log-period Y defined in Eq.(9). Here and throughout,log denotes the natural logarithm (basee), and back- transformations useexp(·). With this convention, Lgeo = exp(Y ) is the geometric mean attractor period for that graph. Differences on theY-scale are multiplicative onLgeo: a change ∆Y corre- sponds to a factor exp(∆Y ) on the geometric mean period [27]. For completeness we also define the arithmetic mean period ¯L in Eq. (8), but the statistical analyses and figures are based onY and Lgeo unless otherwise noted. In the implementation, Y is computed by taking the natural logarithm of each run-specific period Lr, averaging over the M = 100 runs per graph, and storing the

as AvgLogPeriod in the aggregated re- sults file used for the regression. The analysis dataset thus contains one row per graph: the response Y and the associated structural descrip- tors (N, ¯k, MSP, ABC,C ). In a small number of descriptive plots of the raw period distribution we use a base-10 logarithmic axis forLr purely for visualisation (to spread out the long right tail); this display choice does not change the values of Y or affect any model estimates. 9 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint 4.8 Statistical modeling and analysis We summarize how the mean period varies with clustering using a central 80% (10th–90th per- centile) contrast onC: from C0.10 to C0.90. On the natural-log scale this maps directly to a ratio on the period scale [ 27]. The regression is fitted to all graphs; for theC0.10–C0.90contrast we use the fitted clustering slopeβC and the observed span (C0.90−C0.10), define ∆Y = βC (C0.90−C0.10), and report exp(∆Y ) as the adjusted ratio on the period scale. In cross-validation we compute the same contrast separately in each fold using the fold-specificβC. Plots display the fitted curve acrossthe fullobservedsupportof C and highlight the central 80% window used for the contrast. Main regression for mean period.We model the average log-periodY (Eq. (9)) with a Gaussian additive model Y =β0 +βCC +βNN +βk ¯k +βMSP MSP + ti(N, ¯k) + ti(C, ¯k) +ε, (16) where N is the number of nodes, ¯k the aver- age degree, and MSP the mean directed short- est path (Eq. (3)). We model the size-degree and clustering-degree interactions with tensor- product smooths ti(N, ¯k) and ti(C, ¯k) using cubic regression-spline marginals. The basis dimensions forN, ¯k, C and MSP are chosen adaptively based on the number of distinct observed values in each covariate and capped at moderate values (typically between 3 and 10 basis functions per margin) to allow mild nonlinearity without overfitting. Models are fitted by penalised least squares via restricted maximum likelihood(REML)withshrinkageselection, using the mgcv package in R with a Gaussian error distribution. AdjustedC0.10-C0.90contrasts and confidence intervals are computed fromβC as described above and reported on the period scale as exp(∆Y ). Robustness models. As robustness checks, we consider two alternative specifications. First, we replace the linear clustering term by a smooth s(C) while keeping the same interaction struc- ture, to test whether a more flexible, nonlinear relationship between clustering and the log-mean period substantially improves the fit. Second, we fit a more flexible generalized additive model with separate smooth terms forN, ¯k, C, MSP and ABC and the same tensor-product interac- tions, treating this as a fully smooth benchmark. These variants are compared to the main model using both in-sample criteria and out-of-sample performance. To assess the importance of the interaction structure, we also compare the main model to reduced additive models without tensor- product terms, and to models without the ti(C, ¯k) interaction, using likelihood-ratioχ2 tests (via anova.gam in mgcv) and report the correspond- ing test statistics andp-values for the interaction effects. Unadjusted period summaries by clustering. Because the analysis dataset contains a single row per graph (the average natural log-periodY over M=100 runs; Sec. 4.7), we describe how period distributions change with clustering by grouping graphs into deciles ofC and summarising the geometric mean period Lgeo = exp(Y ) within each decile. In particular, we report for each decile the number of graphs, the median, and the interquartile range of Lgeo (Table 3). This gives a simple, nonparametric description of how typical periods and their spread vary across the observed range ofC, without making extra assumptions about the shape of the distribution. Validation and diagnostics. Model adequacy is assessed by 10-fold cross-validation at the graph level (folds contain disjoint sets of graphs), predicted-versus-observed comparisons, calibra- tion by prediction percentiles, residual-versus- fitted plots, and normal Q–Q plots. In cross- validationwefitthemainmodelwithlinear C, the smooth-C variant, and the fully smooth bench- mark on the training folds and compare their out- of-sample performance using root mean squared error andR2 on the held-out graphs (graphs not used to fit the model). To check robustness to modeling choices, we use likelihood-ratioχ2 tests 10 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint to compare the main specification to the smooth- C variant and to reduced models without interac- tion terms; these tests quantify the evidence for nonlinearity inC and for the significance of the interaction structure. 5 Results We study how clustering coefficient variation in small-world graphs shapes the long-run dynamics of synchronous, signed threshold Boolean net- works. Using the Gaussian additive model in Sec. 4 (Eq.(16)), we estimate the effect of clus- tering on the average log-periodY while adjust- ing forN, ¯k, and mean directed shortest path (MSP), and allowing tensor product interactions ti(N, ¯k) and ti(C, ¯k) to capture size-degree and clustering-degree structure. Because changing¯k can drive both clustering and MSP in the Watts- Strogatz model, we evaluate predictors jointly to address collinearity. Raw Spearman correla- tions between Y and each structural descriptor are all positive, ranging fromρ≈0.48 for N and ρ≈0.69 for MSP toρ≈0.94 for mean degree, with ρ≈0.81 for C and ρ≈0.85 for average be- tweenness. These strong positive correlations oc- cur because larger, more highly connected graphs also tend to have higher clustering and shorter paths. For this reason we analyse all predictors together in one model, so that the effect of cluster- ing is interpreted as apartial effect after adjusting for N, ¯k and MSP. Average betweenness centrality (ABC) had only a very small additional partial effect in exploratory fits onceN, ¯k, C, and MSP were included, so ABC is not retained in the main regression. Importantly, the positive raw associa- tion betweenC and Y reverses after adjustment: once N, ¯k, and MSP are held fixed, higher clus- tering is associated withshorter average periods. The remainder of this section first quantifies this adjusted clustering effect, then describes the roles of the other graph properties, model fit and diag- nostics, and finally the unadjusted distributions of periods across clustering deciles. 5.1 Primary effect of clustering on attractor period In the main model the regression term for clus- tering C is linear and negative, and we sum- marise this effect by the coefficientˆβC. Numeri- cally, the fitted coefficient isˆβC =−1.4126 with standard error 0.6988, 95% confidence interval [−2.78,−0.04], andp≈0.044. A 0.10 increase in C multiplies the expected geometric mean at- tractor periodLgeo = exp(Y ) by exp(0.10ˆβC)≈0.87, which corresponds to about a13% reduction on the period scale. We also summarise the effect over the central 80% of the observedC values, from the 10th percentileC0.10 = 0.0000 to the 90th percentileC0.90 = 0.4599. Here C0.10 and C0.90denote the 10th and 90th percentiles of the clustering coefficient across all graphs, so the contrast C0.10→C0.90 compares a typical low- clustering graph to a typical high-clustering graph. Over this range the adjusted ratio on the period scale is exp ( (C0.90−C0.10) ˆβC ) ≈0.52, so moving fromC0.10 to C0.90 shortens the ex- pected geometric mean period by roughly48%. Grouped 10-fold cross-validation at the graph level gives very similar effect sizes: across folds the mean clustering slope is¯βC =−1.65, so a0.10 in- crease inC multiplies the geometric mean period by exp(0.10 ¯βC) = 0.85 with 95% CI[0.76, 0.91], and the interdecile ratio is≈0.49 with 95% CI [0.28, 0.65] (Table 2). Fold-specific slopes are negative in every split, ranging fromβC≈−0.93 to βC ≈−2.79, with the corresponding ratios exp(0.10βC) lying between about0.76 and 0.91 and the interdecile ratios between about0.28 and 0.65, so the shortening effect is present across all foldsratherthandrivenbyafewparticulargraphs. Taken together, the model and cross-validated slopes correspond to about a13-15% shortening in the geometric mean period per+0.10 increase in C. 11 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint The shape of the adjusted clustering effect is illustrated in Fig. 4, which plots the effect ofC on predicted log-period after adjusting forN, ¯k and MSP together with partial residuals forC. Each point shows how far a graph’s observedY lies above or below the fitted model once other covariates are accounted for, while the curves show the estimated conditional effect of clustering from the linear- C model and from a more flexible smooth-C model. Both curves decline steadily with C, and the residual cloud follows the same downward trend, confirming that, conditional on the other graph properties, higher clustering is associated with shorter average log-periods. Addi- tional diagnostics based on individual conditional expectation, accumulated local effects, and the numerical derivative of the smooths(C) term show only mild curvature and a nearly constant negative slope over the observedC range; these robustness checks are described in Sec. 5.5. 5.2 Effects of graph properties While the focus is clustering, the main Gaussian additive model (Eq.(16)) also quantifies how size N, mean degree ¯k, and mean directed shortest path (MSP; Eq.(3)) relate to periods when con- sidered together. In this specification we include linear terms forN, ¯k,C, and MSP plus the tensor- product interactionsti(N, ¯k) and ti(C, ¯k). Table 1 reports the parametric coefficients for the linear terms. In the linear plus tensor model, N, ¯k, and MSP all have positive associations withY (av- erage log-period), whereas clustering has a neg- ative association consistent with Sec. 5.1. The fitted coefficients are ˆβN = 0.0173, ˆβ¯k = 0.7503, ˆβC = −1.413, and ˆβMSP = 0.0555 (Table 1). On the geometric-mean period scale, a +10 increase in N multiplies the expected period by exp(10 ˆβN) ≈1.19, about a 19% increase, a +1 increase in mean degree multiplies it by exp( ˆβ¯k) ≈2.1, more than doubling the typi- cal period, and a +1 increase in MSP multi- plies it by exp( ˆβMSP) ≈1.06, a modest effect compared with degree. For clustering, the per unit ratio exp( ˆβC) ≈0.24 corresponds to the smaller +0.10 contrast discussed above, where exp(0.10 ˆβC)≈0.87, about a13% reduction. As before, we do not compare raw per unit magni- tudes across predictors because units differ; the proportional change column in Table 1 provides a compact way to read off effects on the period scale. The interactionti(N, ¯k) indicates that the in- fluence of size depends on degree: when¯k is near 2 the effect of increasingN is negligible, but at moderate to high¯k the effect ofN becomes clearly positive, consistent with a positive size-degree in- teraction. The ti(C, ¯k) term captures that the influence of clustering depends on degree: at low ¯k changingC has little impact, whereas at higher ¯k increasing clustering pulls the fitted surface to- ward shorter log-periods, reinforcing the negative conditional effect ofC found in Sec. 5.1. These patterns are visualised in Fig. 1, which shows the predicted log(mean period) over( C, ¯k) at fixedN and MSP, and in Fig. 2, which shows the corre- sponding surface over(N, ¯k) with the simulated design points overlaid. Dense, highly connected networks can support long attractor periods, but in this model adding triangles pulls the dynamics back toward shorter cycles. Partial-dependence profilesandtensor-productsurfacesfromthemain linear-plus-tensor model show the same pattern: mean degree has the strongest effect, MSP andN have positive but milder effects, and once these are accounted for the additional contribution of average betweenness is small. To quantify the additional fit provided by the interaction terms, we also compare the main model to reduced ad- ditive models without tensor-product terms and to a model without theti(C, ¯k) interaction using likelihood-ratio χ2 tests (anova.gam in R); the corresponding test statistics andp-values are re- ported in Appendix Table 4, which summarises the statistical significance of the interaction struc- ture. 12 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint Table 1: Parametric coefficients for the mean log-period model.Outcome Y is the average log-period as in Eq. (9); predictors:N, ¯k, C, and MSP. Entries are adjusted per unit effects on the log scale with Wald 95% CIs; the last column reports the proportional changeexp( ˆβ) per +1 unit. ForC,

interpret the prespecified+0.10contrast via exp(0.10ˆβC)≈0.87. Term Estimate Std. Error 95% CI low 95% CI high Ratio per +1 Intercept −2.754 0.150 −3.049 −2.459 - N 0.0173 0.00188 0.0136 0.0210 1.017 ¯k 0.750 0.0552 0.642 0.859 2.12 C −1.413 0.699 −2.78 −0.043 0.24 MSP 0.0555 0.0285 −0.0004 0.111 1.06 Figure 1: Predicted log(mean period) over clustering and mean degree.Contour surface of the fitted linear-plus-tensor model over clustering coefficientC (horizontal axis) and mean degree¯k (vertical axis), with other covariates held at their medians. Periods increase strongly with degree, while at any fixed degree increasingC shifts the surface downward, especially at moderate and high ¯k, illustrating that higher clustering shortens the predicted average period in denser networks. 5.3 Model fit The additive Gaussian model forY captures most of the between-graph variation in mean attrac- tor period. For the primary linear plus tensor model (Eq. (16)) the in-sample performance is strong (adjustedR2 = 0.94, deviance explained 94.3%, RMSE on the log scale≈0.51, scale es- timate≈0.27, n = 330). Standard predicted- versus-observed comparisons, calibration by pre- diction percentiles, residual-versus-fitted plots, and normal Q-Q plots show good agreement and Figure 2: Predicted log(mean period) over size and mean degree with design points. Filled contours of the fitted linear-plus-tensor model over network sizeN (horizontal axis) and mean degree ¯k (vertical axis), with the other covariates held at their medians. Black points mark the simulated (N, ¯k) combinations. Predicted periods increase with bothN and ¯k, and the effect of size is strongest at higher mean degree. The design points show that much of the plotted surface is supported by simulated graphs; however, the upper-left region (small N, higher ¯k) and a few boundary areas lack design points and should be viewed as extrapolation beyond the observed design. no major systematic deviations, supporting the adequacy of the Gaussian GAM on the log scale. These goodness-of-fit diagnostics (predicted vs. observed, calibration by deciles, residuals vs. fit- ted values, and a normal Q-Q plot) are shown in Appendix Figs. 5-8; together they indicate a ro- bust fit with good calibration and approximately Gaussian residuals. When graphs are grouped into deciles of the prediction, mean observed and 13 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint mean predicted log-period agree closely, with dif- ferences generally within about0.07 log units and at most≈0.16, indicating that calibration holds well across the spectrum of fitted values. Grouped 10-fold cross-validation at the graph level shows that out-of-sample performance is similar to the in-sample fit and that the estimated clustering effect is stable across folds. Table 2 summarises the cross-validated slopesβC and the implied period ratios for a+0.10 increase inC and for the interdecile contrastC0.10→C0.90. 5.4 Attractor period distributions by cluster- ing As a descriptive complement to the regression model, we also look directly at the raw data. In this subsection we study theunadjusted distribu- tions of the graph-level geometric mean attractor period Lgeo = exp(Y ) across different values of the clustering coefficientC. We sort graphs byC and split them into ten or- dered groups labelledD01-D10 using a rank-based decile cut onC (graphs with exactly the same C value stay together). Because many graphs haveC≈0, only eight of these deciles contain graphs: D02 is the first non-empty bin (n = 112), and the remaining non-empty bins areD04-D10, which mostly contain about 33 graphs each, ex- cept D04 which has 21 (counts and summaries are given in Table 3). Figure 3 shows violin and box plots ofLgeo by clustering decile, using a log10 y-axis for display. As clustering increases, both the typical period and the spread grow. In the lowest non-empty decile (D02) almost all graphs have very short peri- ods (median = 1.00, interquartile range= 0.000). In the highest decile (D10) the median rises to about 37.3 and the interquartile range is over60 time steps, with a long right tail. The intermedi- atedecilesmovesmoothlybetweentheseextremes. Table 3 summarises the same pattern numerically: median Lgeo rises from 1.00 to≈37.3, and the interquartile range grows from essentially zero to more than 60 time steps. These raw patterns are descriptive, not causal. They show that,without adjustment, graphs with higher clustering tend to have longer and more variable periods, consistent with the positive Spearman correlation betweenC and Y reported earlier. In the Watts-Strogatz model, however, clustering does not change in isolation: higher C is typically accompanied by changes in net- work size, degree, and path structure. The ad- justed analysis that conditions onN, ¯k, and MSP (Sec. 4) shows that thepartial effect of clustering is in the opposite direction: when these other properties are held fixed, higher C shortens the expected geometric mean period (Sec. 5.1). Figure 3: Geometric mean period by clustering decile (log10 y-axis). Violin and box plots of the geometric mean attractor periodLgeo for each clustering decile. As clusteringC increases, both the median and the spread ofLgeo increase, so long typical periods are much more common in the high-clustering groups. The log10 y-axis is used only to spread out the long right tail for visual clarity and does not affect the underlying data or the regression model. 5.5 Robustness We ask whether the estimated clustering effect could be an artefact of our modelling choices, in particular the interaction terms and the assump- tion of a linear effect inC. First, we compare the main linear-plus-tensor model to simpler models using likelihood-ratio χ2 tests (Table 4). Adding the tensor-product terms ti(N, ¯k) and ti(C, ¯k) to an additive model with only linear effects of N, ¯k, C and MSP reduces the deviance by about125.5 units for15.9 14 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint Table 2: Cross-validated clustering effect.Grouped 10-fold cross-validation at the graph level for the main linear-plus-tensor model. For each fold we refit the model, record the clustering slopeβC, and derive the corresponding period ratios for a+0.10increase inC and for the interdecile contrastC0.10→C0.90. The table reports the mean across folds, the standard deviation (SD), and an empirical 95% interval (2.5th-97.5th percentiles). Quantity Mean SD 2.5% 97.5% βC −1.65 0.52−2.61−0.79 Ratio for +0.10in C 0.85 0.04 0.76 0.91 Ratio forC0.10→C0.90 0.49 0.10 0.28 0.65 Table 3: Decile sizes and summaries of geometric mean period by clustering decile.n is the number of graphs; medians and interquartile ranges (IQRs) are forLgeo; Cmid is the bin midpoint (3 s.f.). Decile n Median Lgeo IQR Lgeo Cmid D02 112 1.00 0.000 0.000 D04 21 1.07 0.466 0.0239 D05 32 2.54 2.92 0.0485 D06 33 6.55 29.3 0.0795 D07 33 6.91 60.0 0.130 D08 33 7.62 26.5 0.241 D09 33 16.1 22.9 0.388 D10 33 37.3 66.2 0.530 effective degrees of freedom (p<0.001). Adding the ti(C, ¯k) interaction on top of a model that already includes ti(N, ¯k) reduces deviance by a further 68.0 units for 12.8 degrees of freedom (p < 0.001). These tests support keeping both interactions and show that the effect of clustering depends on degree in a statistically important way. Second, we then test whether a more flexible, nonlinear term in clustering improves the fit. We replace the linearC term by a smooth function s(C) in an otherwise identical Gaussian additive model (Sec. 4, Eq.(16)), keeping the same linear terms forN, ¯k and MSP and the same tensor- product interactions ti(N, ¯k) and ti(C, ¯k). The smooth s(C) is a cubic spline with basis sizekC, and it is penalised so that the fitted curve is kept as smooth as possible and does not introduce un- necessary wiggles. The fitted effective degrees of freedom fors(C) are only slightly above1, which means the estimated effect ofC is almost linear, with only mild curvature. Deviance explained and adjustedR2 are essentially unchanged com- pared with the linear-C model, and the smooth-C variant has a slightlyhigher residual deviance (Ta- ble 5; ∆ deviance≈−0.14 at ∆df≈0.29). There is therefore no evidence that a nonlinear effect in C improves the fit, and we retain the linear specification forC. The shape of the fitted clustering effect is also consistent with this conclusion. Figure 4 plots partial residuals forC from the linear model to- gether with the estimated contributions of the linear C term and the smooths(C) term. Both fitted curves decline steadily withC, they are almost straight over the observed range, and the residual cloud follows the same downward trend. This indicates an approximately constant nega- tive slope on the log scale as clustering increases, and shows that the linear term inC is a parsimo- nious and adequate summary of the conditional effect. Finally, we fitted a more flexible comparison model with separate smooth terms forN, ¯k, C, MSP and ABC, plus the same tensor-product interactions. This fully smooth model gives simi- lar deviance explained and out-of-sample perfor- mance to the main linear-plus-tensor specification. Partial-dependence profiles for the additional co- variates (Appendix Fig. 9) show that mean de- gree has the strongest effect onY, MSP andN 15 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint have positive but milder effects, the partial effect of ABC is small once the other descriptors are included, and the negative effect ofC remains, consistent with the main model. Figure 4: Effect of clustering after adjustment and partial residuals.Partial residuals forC from the linear model (grey points), plotted againstC, together with the estimated conditional effect of clustering from the linear-C model (green curve) and from the smooth-C model (blue curve), each with a 95% confidence band. The residual cloud follows a clear downward trend, and both fitted curves decline steadily withC, confirming a negative conditional association between clustering and the average log-period after adjusting forN, ¯k and MSP. Reproducibility. All code and data will be pub- licly available. A tagged GitHub release will in- clude scripts to rebuild all tables and figures, the processed data, a manifest (commit hash and ran- dom seeds), and an environment file with package versions. We will archive the release with a DOI and cite it; until the DOI is issued, materials are available on request. 6 Discussion We quantify how small-world wiring controls cycle length in synchronous, signed threshold Boolean networks. Adjusting for size N, mean degree ¯k, and mean directed shortest path (MSP) in a single model with size-degree and clustering- degree tensor terms (Sec. 4), higher clustering C is linked to shorter attractor periods. In the main regression, a0.10 step inC multiplies the expected geometric mean periodLgeo = exp(Y ) by exp(0.10 ˆβC)≈0.87 (about 13% shorter), and moving from C0.10 = 0.0000 to C0.90 = 0.4599 yields an average≈48% reduction (exp(∆Y )≈ 0.52; Sec. 5). Cross-validated estimates are very similar: grouped 10-fold CV at the graph level gives a mean slope¯βC≈−1.65, so a0.10 increase in C multiplies the geometric mean period by exp(0.10 ¯βC)≈0.85 (Table 2). The shortening effect of clustering therefore generalises across held-out graphs and corresponds to about a13- 15% shortening in the geometric mean period per +0.10increase inC. Because the WS rewiring probability alters tri- angle density and path length together,C and MSP co-vary. We therefore identify the effect of clustering by modelling therealised graph statis- tics and estimating thepartial association ofC with the average log-periodY while conditioning on MSP (and onN and ¯k). The negative slope for C persists under this joint adjustment, under cross-validation, and when the linearC term is re- placed by a penalised smooths(C). The smooth adds only slightly more than one effective degree of freedom and does not materially change the fitted effect or improve fit in terms of deviance explained, R2, or information criteria (Appendix Table 5). In Fig. 4 the fitted linear and smooth curves are almost indistinguishable, reinforcing that a linear term inC is adequate. We therefore interpret ˆβC as a robustconditional association between clustering and shorter periods within the WS model. Mechanistically, the pattern aligns with the ordered-chaotic spectrum for Boolean dynamics [3, 9, 2, 28]. Triangles create short feedback loops and local redundancy; together these features lower effective sensitivity, slow damage spreading, and shorten transients [25, 26, 35]. Empirically, the fitted curves forC in Fig. 4 are close to linear and remain negative over the observed range, sug- gesting that within our WS range clustering acts as a steady, order-promoting factor rather than a sharp switch. This complements prior work 16 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint on information storage and perturbation spread- ing by providing a calibrated effect size on cycle length while controlling for size, degree, and path structure. Other structural levers behave coherently. Holdingtheremainingpredictorsattheirmedians, increasing ¯k by +1 multiplies the expected geo- metric mean period byexp(0.750)≈2.12, adding +10 nodes multiplies it by exp(10×0.0173)≈ 1.19, and MSP is positively associated (ratio exp(0.0555)≈1.06 per unit) (Table 1). The fitted interaction betweenN and ¯k indicates a positive size-degree interaction: increases in N matter most at moderate to high¯k. The corresponding clustering-degree interaction shows that the in- fluence of clustering depends on degree: at low ¯k varyingC has little impact, whereas at higher ¯k increasing C moves predictions toward shorter log-periods, consistent with a negative conditional effect of clustering. Our study has a specific scope: we work with directed, signed WS graphs with synchronous up- dates and strict threshold rules over the ranges reported in Sec. 4. Model fit and calibration are strong in grouped 10-fold CV at the graph level (Sec. 5.3, Table 2), but extrapolations to asyn- chronous updates [29], other logic families (e.g. canalising or nested-canalising) [21, 37], heavy- tailed degree sequences [32, 14, 6], or strong mod- ularity [33] should be made with care. In sum, within the WS model analysed, clus- tering exerts a clear, independently estimable, monotone shortening effect on attractor periods once size, degree, and path length are accounted for. This links a concrete structural motif, tri- angles, to a quantitative, predictive handle on long-run dynamics in threshold networks relevant to synthetic circuits and interpretable models of cellular decision making [19, 5, 38]. 7 Conclusion Clustering is a simple, interpretable wiring fea- ture with a clear dynamical consequence in signed, threshold Boolean systems. In Watts-Strogatz small-world graphs with synchronous updates, the partial effect of global clustering on the av- erage log-periodY is monotone and close to lin- ear, remains negative over the observed support (0≤C≤0.6274), and persists after conditioning on N, ¯k, and MSP. These findings are consis- tent with theory on information storage, damage spreading, and loop closure. Because the effect is expressed in units ofC itself, it is portable to empirical networks where clustering is measured directly. Practically, clustering joins degree and size, as well as path-length structure captured by MSP, as a quantitative control for tuning conver- gence speed and the length of repeating cycles in discrete dynamical systems.

Maram Alqarni is supported by a scholarship fromtheSaudiArabianCulturalMission(SACM). This research was supported by the Australian Re- search Council through the Australian Research Council Centre of Excellence for Plant Success in Nature and Agriculture (CE200100015).

[1] Carlos Gershenson. Introduction to ran- dom boolean networks. arXiv preprint nlin/0408006, 2004. [2] Barbara Drossel. Random boolean networks. Reviews of nonlinear dynamics and complex- ity, pages 69–110, 2008. [3] S. A. Kauffman. Metabolic stability and epigenesis in randomly constructed genetic nets. Journal of Theoretical Biology, 22(3): 437–467, 1969. [4] Martin Hopfensitz, Christoph Müssel, Markus Maucher, and Hans A Kestler. At- tractors in boolean networks: a tutorial. Computational Statistics, 28:19–36, 2013. [5] Assieh Saadatpour and Réka Albert. Boolean modeling of biological regulatory networks: a methodology tutorial.Methods, 62(1):3–12, 2013. 17 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint [6] Maximino Aldana. Boolean dynamics of net- works with scale-free topology.Physica D: Nonlinear Phenomena, 185(1):45–66, 2003. [7] Konstantin Klemm and Stefan Bornholdt. Stable and unstable attractors in boolean networks. Physical Review E, 72(5):055101, 2005. [8] Phillip Bonacich. Some unique properties of eigenvector centrality. Social networks, 29 (4):555–564, 2007. [9] Ilya Shmulevich and Stuart A Kauffman. Ac- tivities and sensitivities in boolean network models. Physical review letters, 93(4):048701, 2004. [10] Stefan Bornholdt and Kim Sneppen. Robust- ness as an evolutionary principle.Proceed- ings of the Royal Society of London. Series B: Biological Sciences, 267(1459):2281–2286, 2000. [11] Fangting Li, Tao Long, Ying Lu, Qi Ouyang, and Chao Tang. The yeast cell-cycle net- work is robustly designed. Proceedings of the National Academy of Sciences, 101(14): 4781–4786, 2004. [12] Reinhard Diestel. Graph Theory. Springer, 5th edition, 2017. [13] Frank Harary. Graph Theory. Addison- Wesley, 1969. [14] Réka Albert and Albert-László Barabási. Statistical mechanics of complex networks. Reviews of modern physics, 74(1):47, 2002. [15] Mark EJ Newman. The structure and func- tion of complex networks.SIAM review, 45 (2):167–256, 2003. [16] Mikaela Koutrouli, Evangelos Karatzas, David Paez-Espino, and Georgios A Pavlopoulos. A guide to conquer the biologi- cal network era using graph theory.Frontiers in bioengineering and biotechnology, 8:34, 2020. [17] Duncan J Watts and Steven H Strogatz. Col- lective dynamics of ‘small-world’ networks. Nature, 393(6684):440–442, 1998. [18] Mark EJ Newman. The mathematics of net- works. The new palgrave encyclopedia of economics, 2(2008):1–12, 2008. [19] Leon Glass and Stuart A Kauffman. The logical analysis of continuous, non-linear bio- chemical control networks.Journal of theo- retical Biology, 39(1):103–129, 1973. [20] S. Karanam and W.-J. Rappel. Boolean mod- elling in plant biology.Journal of Theoretical Biology, 2022. [21] Stuart Kauffman, Carsten Peterson, Björn Samuelsson, and Carl Troein. Genetic net- works with canalyzing boolean rules are al- ways stable. Proceedings of the National Academy of Sciences, 101(49):17102–17107, 2004. [22] Tiago P Peixoto and Barbara Drossel. Boolean networks with reliable dynamics. Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, 80(5):056102, 2009. [23] Hadeel Kittaneh, Filippo Castiglione, and Abdul Salam Jarrah. Stability of threshold boolean networks. Journal of Complex Net- works, 13(4):cnaf011, 2025. [24] Xiao Yang, Nilam Ram, Peter CM Mole- naar, and Pamela M Cole. Describing and controlling multivariate nonlinear dynamics: A boolean network approach.Multivariate Behavioral Research, 57(5):804–824, 2022. [25] Joseph T Lizier, Siddharth Pritam, and Mikhail Prokopenko. Information dynamics in small-world boolean networks.Artificial life, 17(4):293–314, 2011. 18 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint [26] Qiming Lu and Christof Teuscher. Damage spreading in spatial and small-world random boolean networks. Physical Review E, 89(2): 022806, 2014. [27] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.An introduc- tion to statistical learning: with applications in R, volume 103. Springer, 2013. [28] Florian Greil and Kevin E Bassler. Attrac- tor period distribution for critical boolean networks. Physical Review Letters, 102(4): 048101, 2009. [29] Albert I. Albert R. Saadatpour, A. Attractor analysis of asynchronous boolean models of signal transduction networks. Journal of Theoretical Biology, 266(4):641–656, 2010. [30] Deok-Sun Lee. Evolution of regulatory net- works towards adaptability and stability in a changing environment.Physical Review E, 90(5):052822, 2014. [31] Jason Lloyd-Price, Abhishekh Gupta, and Andre S Ribeiro. Robustness and informa- tion propagation in attractors of random boolean networks. PloS one, 7:1–6, 2012. [32] Albert-László Barabási and Réka Albert. Emergence of scaling in random networks. Science, 286(5439):509–512, 1999. [33] Neeraj Pradhan, Subinay Dasgupta, and Sitabhra Sinha. Modular organization en- hances the robustness of attractor network dynamics. Europhysics Letters, 94(3):38004, 2011. [34] Caitlyn Parmelee, Samantha Moore, Kather- ine Morrison, and Carina Curto. Core motifs predict dynamic attractors in combinatorial threshold-linear networks.PloS one, 17(3): e0264456, 2022. [35] Carlos Handrey A Ferraz and Hans J Her- rmann. The kauffman model on small-world topology. Physica A: Statistical Mechanics and its Applications, 373:770–776, 2007. [36] Ugo Bastolla and Giorgio Parisi. Closing probabilities in the kauffman model: An an- nealed computation. Physica D: Nonlinear Phenomena, 98(1):1–25, 1996. [37] Maximino Aldana and Philippe Cluzel. A natural class of robust networks.Proceedings of the National Academy of Sciences, 100 (15):8710–8714, 2003. [38] Jorge GT Zañudo and Réka Albert. Cell fate reprogramming by control of intracellu- lar network dynamics.PLoS computational biology, 11(4):e1004193, 2015. [39] Andy Beatty, Christopher R Winkler, Thomas Hagen, and Mark Cooper. Predict- ing trait phenotypes from knowledge of the topology of gene networks. bioRxiv, page 2021.06.29.450449, 2021. [40] Linton C Freeman. A set of measures of centrality based on betweenness.Sociometry, pages 35–41, 1977. [41] Andrew Pomerance, Edward Ott, Michelle Girvan, and Wolfgang Losert. The effect of network topology on the stability of discrete state models of genetic control.Proceedings of the National Academy of Sciences, 106 (20):8209–8214, 2009. 19 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint Appendix A Model comparisons and robustness A.1 Likelihood-ratio tests for interaction terms Here we summarise the likelihood-ratio (LR)χ2 tests that compare the main linear-plus-tensor model with simpler models. The LR statistic is the change in deviance (χ2) with approximate de- grees of freedom (df);p-values are from likelihood- ratio χ2 tests (using anova.gam in R). The first comparison shows that adding both interaction terms ti(N, ¯k) and ti(C, ¯k) gives a large reduction in deviance (χ2≈125.5 for about 15.9 degrees of freedom,p < 0.001) compared with an additive model with only linear effects of N, ¯k,C, and MSP. The second comparison shows that adding ti(C, ¯k) on top ofti(N, ¯k) alone also gives a strong improvement in fit (χ2≈68.0 for about 12.8 degrees of freedom,p< 0.001). These tests support the use of the interaction structure in the main model. A.2 Linear vs smooth clustering effect We also compared the main model with a linear term in the clustering coefficientC to a variant where C is replaced by a penalised splines(C). The smooth-C model has a slightly larger residual deviance and almost the same degrees of freedom, so there is no evidence that allowing extra curva- ture inC improves the fit. The change in deviance is negative (∆ deviance ≈ −0.14 for a small change in df), so the smooth-C model does not give a better fit. This supports keeping a simple linear term inC in the main specification. A.3 Model diagnostics In this section we show the main goodness-of-fit diagnostics for the Gaussian additive model for Y (average log-period). All logarithms here are natural logs (log basee). A.4 Predicted vs observed and calibration Figure 5 shows the predicted versus observed log- period for all graphs, with a simple linear trend line. Points lie close to the diagonal, indicating a good overall fit. Figure 6 shows calibration by prediction deciles: graphs are grouped into ten bins by their predicted value, and we plot mean observed versus mean predicted log-period in each bin. Figure 5: Predicted vs observed average log-period. Scatter plot of predicted versus observedY = log(period) for all graphs under the main linear-plus-tensor model, with a fitted straight line. Points lie close to the diagonal, indicating good overall fit. Figure 6: Calibration by prediction deciles. Mean observed versus mean predicted log-period by deciles of the predicted value. Point size shows the number of graphs in each decile. Most points lie close to the diagonal, and differences are small across the range, indicating good calibration. A.5 Residual diagnostics Figure7plotsdevianceresidualsagainstthelinear predictor. There is no strong pattern and the smooth trend line is close to zero across the range, suggesting that the mean structure is adequate. Figure 8 shows a normal Q-Q plot of the residuals. Points stay close to the reference line, except for 20 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint Table 4: Likelihood-ratio tests for interaction terms.Each row compares a reduced model with the main linear-plus-tensor model forY (Eq. (16)). The first comparison tests the joint contribution of both interaction termsti(N, ¯k) and ti(C, ¯k). The second comparison tests the additional contribution of the clustering-degree interactionti(C, ¯k) when ti(N, ¯k) is already included. The LR statistic is the change in deviance (χ2) with approximate degrees of freedom (df). Comparison Full model Reduced model Terms tested (dropped) χ2 df p-value Additive vs. full Main model: linear terms + ti(N, ¯k) + ti(C, ¯k) Additive model: linear terms only (no ti terms) ti(N, ¯k), ti(C, ¯k) 125.5 15.9 < 0.001 No ti(C, ¯k) vs. full Main model: linear terms + ti(N, ¯k) + ti(C, ¯k) Main model without ti(C, ¯k) ti(C, ¯k) 68.0 12.8 < 0.001 Table 5: Linear vs smooth clustering effect.Comparison of the main model with a linear term in clustering coefficientC and a variant with a penalised splines(C). The smooth-C variant has slightly higher residual deviance, so the linear term is preferred. Model Residual df Residual deviance ∆df ∆deviance Linear C 309.09 84.21 - - Smooth s(C) 308.80 84.36 0.29 −0.14 mild deviations in the far tails, so the Gaussian assumption on the log scale is reasonable. Figure 7: Residuals vs linear predictor. Deviance residuals of the main model versus the linear predictor. The smooth curve is close to zero and there is no strong structure, indicating that the model captures the main trends inY. Figure 8: Normal Q-Q plot of deviance residuals. Normal Q-Q plot for the deviance residuals of the main model. Points are close to the

line, with mild deviations in the tails, consistent with (though not definitive evidence for) a Gaussian error model on the log scale. A.6 Additional partial-dependence plot Figure 9 shows partial-dependence profiles for N, ¯k, MSP, and average betweenness centrality 21 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint (ABC), holding the other covariates at their me- dian values. These plots illustrate that mean degree has the strongest effect onY, MSP and N have positive but milder effects, and the ad- ditional contribution of ABC is small once the other variables are included. Figure 9: Partial-dependence profiles for graph descriptors. Partial-dependence curves for Y = log(period) as a function of each covariate (N, ¯k, MSP, and ABC), with other covariates fixed at their medians. Shaded bands are 95% confidence intervals. Mean degree¯k has the strongest effect, MSP andN are positive but milder, and the effect of ABC is small once the other descriptors are included. 22 .CC-BY 4.0 International licenseavailable under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprintthis version posted January 2, 2026. ; https://doi.org/10.64898/2026.01.01.697270doi: bioRxiv preprint