Papers

Continuous-time Markov chains (CTMCs) provide the backbone for modeling discrete stochastic dynamics across applied, physical, and biological sciences. Their integration with modern gradient-based machine learning, however, is limited by the hard categorical event selection intrinsic to Gillespie-type simulation algorithms. We exploit the affine state update to obtain the exact one-step conditional-mean sensitivity by differentiating normalized reaction propensities. We pair this backward rule with exact forward trajectories to define the propensity straight-through (PST) estimator. At the trajectory level, we show that one-step sensitivities composed across events can depart from the exact multistep sensitivity. We derive the resulting per-step discrepancy in closed form and prove that it vanishes identically for affine downstream dependence. PST matches the accuracy of Gumbel-Softmax straight-through across all benchmarks: reversible dimerization (0.06% error), a genetic oscillator (1.7% error), a 50-task repressilator suite (0.17% median error), and patch-clamp ion-channel recordings ($R^2$ = 0.988). Under matched settings, PST converges 3.0-fold faster on the oscillator and 2.1-fold faster on the ion channel. At deep-learning scale, PST trains a 203,796-parameter stochastic reaction network with hard sampling, reaching 98.22% MNIST digit classification accuracy. By differentiating an exact conditional mean rather than a relaxed sample, PST offers a temperature- and Gumbel-free path to scalable gradient-based learning through exact stochastic trajectories.

Strongly connected components (SCCs) characterize modular structure in directed networks but are fragile to single node failures. We study strongly biconnected components (SBCs), which are the set of nodes in which every node pair remains mutually reachable after the removal of any single node, as a more robust notion of connectivity. Using a generating function formalism, we derive the size of the giant SBC and analyze its percolation behavior under random node and link removal. We show that the giant SBC emerges at the same threshold as the giant SCC but grows more slowly due to stricter connectivity requirements. We also applied our theoretical framework to real-world biological networks including gene regulatory networks and neural connectome. Our framework provides insight into the interplay between connectivity, redundancy, and robustness in complex directed systems.

Computational pathway design often focuses on evaluating selected pathways or optimizing fluxes in a fixed network, but gives less direct access to the combinatorial question of which other enzyme subsets of the network can support productive alternative pathways. A structured computational analysis of these networks can act as a valuable pre-step to the pathway design process. We present here a systematic approach for exploring biochemical pathway alternatives across enzyme subsets, using a computational methodology based on a rule-based modeling of the enzymes: for a biochemical system with enzyme set $S$, we evaluate all subsets $s \subseteq S$ by generating chemical reaction spaces, searching for integer-hyperflow pathways from prescribed inputs to target products, and organizing feasible subsets by set inclusion. This yields an inclusion-ordered landscape of pathway feasibility and carbon efficiency. We apply the approach to the non-oxidative pentose phosphate pathway, to non-oxidative glycolysis, and to glycolysis. Across these systems, feasible subsets occupy only a moderate fraction of all enzyme subsets, but the structure of this feasible region differs strongly between the systems. Larger enzyme sets do not consistently improve carbon efficiency when every enzyme in the tested subset is required to participate in the pathway. Instead, performance depends on specific enzyme combinations. The resulting subset landscapes are valuable means for identifying essential enzymes, candidate redundancies, and small high-performing enzyme subsets. By making the enzyme-subset landscape itself the object of analysis, the approach addresses the gap between detailed evaluation of individual candidate pathways and early-stage design decisions about which enzyme combinations are worth investigating at all.

Conventional annotation of single-cell RNA-sequencing (scRNA-seq) data relies heavily on manual, marker-based thresholding, an approach that can obscure subtle transcriptomic gradients and collapse functionally distinct cell states into broad, heterogeneous populations. Here we apply the Gaussian multi-Graphical Model (GmGM) framework, which jointly infers cell-cell and gene-gene dependency structure from a single scRNA-seq data matrix, to a 10x Genomics PBMC dataset. Ten independent GMGM-Leiden clustering runs were integrated into a robust consensus partition using a soft cluster ensemble approach and benchmarked against reference cell-type annotations. This strategy yielded stable cluster partitions that resolve biologically meaningful sub-populations not distinguished by the reference annotation. In parallel, for each cluster, gene co-expression modules were extracted from the fitted model via consensus Leiden clustering across resolutions, evaluated using standard network metrics, and validated functionally with the Network Enrichment Analysis Test (NEAT), which confirmed non-random enrichment signal. A module-scoring procedure linked network topology to per-cell, per-cluster expression signatures, and a novel extension of GmGM, recovering a shared cell-cell network together with population-specific gene networks in a single model run, was demonstrated in a case study on the CD4+ T-cell population. These results indicate that GmGM provides a unified, reproducible framework for joint cell clustering and gene-network inference, capable of revealing cellular structure beyond that captured by conventional pipelines.

Reengineered biological systems have the potential to revolutionize chemical and material production, enhance critical mineral recovery, serve as threat sensors and improve human health. Unfortunately, the extreme complexity of organisms has made it difficult to achieve this potential in all but the simplest cases. Recent technological advances, however, have provided a foundation for solving this problem. Here, we describe a capability for accelerating the reengineering of cells by providing accurate predictions of the impact of genetic or environmental changes on cell phenotype. This digital phenome platform (DigiPhen) consists of integrated experimental, analytical and modeling workflows for building a digital representation of microbial or plant systems. It is designed around an expanding set of interchangeable, interconnecting software and experimental modules that can accurately represent the mechanistic determinants of phenotype. The DigiPhen platform will systematically collect data on cell composition, spatial organization, metabolic pathways and regulatory networks in a semi-autonomous fashion and use this information to build modular, multi-scale models of biological systems. These models will be used to predict molecular and environmental changes needed for producing desired biological outcomes. DigiPhen is intended to be the heart of community research campaigns that will meet the immediate needs of individual researchers while fulfilling long-term goals of the scientific community. Altogether, the DigiPhen platform represents a new paradigm for building predictive models of biological systems.

A composite structural index summarises a network in one number, and for a triangle-based index it is spectrally redundant: Tr(A^3) is the third moment of the adjacency spectrum. The non-redundant content sits one level down, in diag(A^3), which depends on eigenvectors and is not spectrally determined. A corollary in the theory paper predicted that the global scalar should tie sharpened spectral baselines rather than beat them, while the node-wise attribution should do better where the number of structural epicentres is unknown. We construct Omega-N by localizing each of the four factors. The direct localization is badly conditioned; two corrections from published practice fix it, a configuration-null excess per factor and a personalized-PageRank neighbourhood at several scales, giving ten interpretable features per node, with no attributes, training or embeddings. Against a recursive feature engine at five levels of recursion, Omega-N wins on three and ties on two of the six in-domain evaluations, the sixth a declared null where every arm returns chance, with ten features against its 28 to 252 before pruning. Two statistics from the graph and labels, not from performance, partition the eight benchmarks without error, and the two they exclude are the two on which it loses. The strongest application is drug-target prioritisation on protein interaction networks: +0.032 to +0.103 AUPRC over a six-feature centrality battery and +0.084 to +0.208 over the four-feature one, across three constructions, replicated on an independent AP-MS network and label source (degree-matched: +0.0723 on STRING, +0.0560 on BioPlex, p=0.00195). Adding Omega-N to centralities plus Node2Vec changes nothing. The claim is narrow and it is the point: ten named features, computed without training, match or beat hand-crafted centralities and a recursive engine, and do not touch learned representations.

Network Hamiltonian Models (NHMs) provide an efficient framework for modeling the aggregation of interacting particles (e.g., the condensation of proteins into gel-like, oligomeric, or fibrillar states), representing the system as a network whose edges represent bound interactions. Terms within the network Hamiltonian represent multi-body interactions governing aggregation behavior, and are specified via topological degrees of freedom. Because motional degrees of freedom are not explicitly represented within the NHM, their influence must be indirectly accounted for by introduction of corresponding terms to the Hamiltonian. Here, we describe specifications for terms representing motional degrees of freedom of two, three, and four-body interactions. We also consider the impact of these terms on the aggregation states of a minimal system governed only by a pairwise edge potential, showing that three and four-body interactions favor condensation of the system into small droplet-like structures at low temperature.

Metabolic networks play pivotal roles in understanding the evolution of organisms, microbiome dynamics and disease prevention and treatment. This study presents a comprehensive analysis of metabolic network properties across 10912 organisms spanning Bacteria, Archaea, and Eukarya domains. A novel method for the network construction is introduced, emphasizing the chemical transformations of metabolites. Unlike conventional approaches that link every substrate to every product in all chemical reactions, this method aligns more closely with metabolic pathways, linking products only to their generating substrates. Emphasis is placed on investigating the community structure within metabolic networks, employing the Surprise quality function that better addresses some limitations of previous approaches. Through comparisons between real and randomized versions of the networks, we identify characteristics that cannot be explained solely by the network degree distribution, thus identifying topological properties and community structures that arise from additional evolutionary pressures. This highlights the need for evolutionary models to incorporate additional mechanisms beyond the replication of the degree distribution.

Boolean networks provide a compact framework for modeling regulatory systems, yet their rapidly expanding state spaces make systematic dynamical analysis challenging. Here, we systematically enumerate all non-isomorphic two-node signed regulatory subgraphs with their admissible Boolean update rules and exhaustively characterize their state-transition graphs (STGs) under both synchronous and asynchronous updating. By extracting quantitative dynamical descriptors, we move beyond static wiring diagrams toward a behavior-based classification of network structures. The synchronous analysis identifies deterministic fixed-point and cyclic attractor regimes, while the asynchronous analysis captures one-node-at-a-time dynamics through terminal strongly connected components of the asynchronous STGs. Comparing the two update schemes distinguishes dynamical behaviors that are robust across update assumptions from those that depend on synchronous updating. Despite the diversity of graph-rule combinations, many realizations converge to a limited repertoire of dynamical classes, revealing substantial redundancy between structural and dynamical representations. We further propose a testable framework for investigating whether such local dynamical classes can be used to characterize selected subnetworks within larger Boolean models. To facilitate exploration and reproducibility, we developed BORNA, an interactive web application for exploring the complete network catalogue, STGs, and synchronous and asynchronous simulations, available at https://jafarilab.github.io/BORNA/. Together, these results provide a systematic reference for minimal Boolean network dynamics and a foundation for extending behavior-based analysis to larger regulatory systems.

A fundamental problem in the algebraic study of biochemical reaction networks is to characterize degeneracy. Two-dimensional zero-one networks, where each reactant appears with stoichiometric coefficient zero or one, constitute the smallest biologically relevant class capable of exhibiting degeneracy. In this paper, we provide a complete classification: a two-dimensional zero-one network without trivial species is degenerate if and only if it is a consistent subnetwork of a species refinement of one of two prototypical networks, namely the complete paired-exchange network or the catalytic-pair conversion network. Equivalently, in matrix-theoretic terms, degeneracy is determined entirely by the row patterns of the stoichiometric and reactant matrices, and can be verified by purely structural inspection without computation. Remarkably, every such degenerate network has a steady-state system consisting entirely of binomials, so its steady-state variety is toric. These results also bear on absolute concentration robustness: since any network exhibiting this property necessarily contains a degenerate subnetwork, the minimal degenerate networks characterized here serve as fundamental building blocks for constructing such networks.

In the transcriptional networks of Escherichia coli and Saccharomyces cerevisiae, the eight feed-forward loop (FFL) motifs occur at markedly different frequencies. Although previous studies have linked the abundant C1- and I1-FFLs to specific dynamical functions, a common account of the broader pattern is lacking. To address this gap, we adopt an information-theoretic approach. An FFL transmits information through two paths that share an input and converge on an output, suggesting that their interaction may contribute to differences among motifs. To investigate this, we decompose input-output mutual information (MI) into pathway MI and interference MI (IMI). Here, the pathway MI collects the contributions of individual pathways, whereas IMI identifies the contribution that arises from their interference. We find that within a representative parameter regime, IMI values of FFLs produce a hierarchy that closely resembles the abundance ordering in E. coli, whereas total and pathway MIs do not. Motivated by this observation, we perform constrained optimization that identifies parameter regimes compatible with abundance-like targets for both S. cerevisiae and E. coli. We relate the optimized IMI hierarchy to pathway interference strength and local pathway sensitivities to understand the underlying biophysical origin. Our framework thus probes the correspondence between pathway interference, information transmission, and motif-abundance hierarchies.

Biochemical reaction networks have become a central theoretical framework for implementing molecular computation. A key challenge is finite time computational accuracy, as computation outputs are encoded in limiting steady states (LSSs) of species concentrations while practical implementations operate for only finite time. This work proposes a concise characterization of convergence rate for biochemical computational modules with multiple output species, and rigorously establishes it as being bounded by the eigenvalue with largest (least negative) real part of the Jacobian matrix. Two numerical examples illustrate how the theoretical lower bound shapes the convergence rate range and reveals its dependence on reaction rate constants. This formulation enables systematic evaluation of biochemical computation speed and provides a practical design measure for constructing high-accuracy and error-controlled biochemical computational modules.

Gene expression is an intrinsically stochastic process that generates phenotypic heterogeneity within genetically identical cell populations. While the exact statistical moments of the protein count can be obtained for a broad range of complex models, the corresponding distributions are significantly harder to obtain and intractable in many cases. The classical three-stage model of gene expression, which predicts fluctuations in protein levels as a function of promoter switching, transcription, translation, and degradation events all occurring with linear propensities, illustrates this perfectly; deriving its exact protein distribution remains elusive. Here, using the partitioning property of time-inhomogeneous Poisson processes, we develop an exact mapping of the three-stage model onto a simplified model. The simplified model allows us to formulate two analytical approximations for the full protein distribution of the three-stage model based on a beta-mixture representation of the exact solution for a simpler model. We show that the two approximations are asymptotically exact in different limiting cases and verify their accuracy against simulations for a broad range of parameters. Although approximate, these are the first analytical expressions for protein distributions for the three-stage model that are highly accurate in intermediate regimes.

Biological recognition rarely rests on one strong bond. It works by forming many weak ones at once, between crowded, deformable surfaces in water. This review develops that process as a problem in statistical mechanics. Counting the ways two multivalent objects can bind proves to be the classical monomer-dimer problem on a graph, with a rigorous consequence: the apparent switching of multivalent binding is always a smooth crossover, never a phase transition. Three constraints follow. A repulsive surface layer is obligatory rather than a design choice; bonds do not act independently; and since free energies enter rates exponentially, small changes in receptor number shift binding lifetimes by orders of magnitude. One set of equations then covers antibodies, lipoproteins, and T cell recognition. In each, what decides the outcome is not the strength of any single bond but how a fixed total is spread over many: affinity is a property of a molecule, selectivity a property of an assembly.

5-aminolevulinic acid (5-ALA) promotes fluorescence-based resection of glioblastoma via protoporphyrin IX (PpIX); however, there is considerable heterogeneity in fluorescence intensity, limiting margin distinction. Cell-line studies show that stimulation of the MAPK pathway results in decreased PpIX levels through increased elimination by ABCB1, an efflux transporter, and ferrochelatase (FECH), the enzyme that converts PpIX to heme. However, the mechanism has yet to be studied in human tissue. Using two independent, purity-adjusted sets of primary IDH-wildtype glioblastoma specimens (TCGA-GBM, n=140, discovery; CGGA, n=87, validation), no reproducible correlation between MAPK signaling and its proposed downstream effectors (ABCB1 and FECH) could be detected. Instead, MAPK activity displayed a reproducibly negative relationship with PPOX, the enzyme that converts protoporphyrinogen IX to protoporphyrin IX: TCGA-GBM (n=133 with purity estimates; rho = -0.40, 95% bootstrap CI [-0.54, -0.25], adjusted p < 0.001), CGGA (rho = -0.28, 95% bootstrap CI [-0.47, -0.07], adjusted p = 0.0342). The finding was robust to permutation testing and to alternative approaches to purity adjustment. Here we report, for the first time, replicated human-tissue evidence for this association. If this correlation reflects causality, MAPK regulation of PpIX would act at the point of synthesis rather than at the clearance steps implied by previous cell-culture studies. We observed differential activation of the MAPK pathway with increased activation in the classical subtype and decreased activation in the proneural subtype. There was no correlation between MAPK pathway activation and overall survival (log-rank p = 0.73). The association between PPOX expression and MAPK pathway activation is biological rather than prognostic, identifying the MAPK pathway as a candidate handle for increasing 5-ALA fluorescence.

Boolean networks are widely used to model gene regulatory Attractors of Boolean networks model stable gene-expression patterns, yet deriving their properties from network structure remains an open problem. We solve this problem for critical $K=1$ networks by showing that their feedback loops induce a permutation whose order bounds the average attractor length both above and below by universal constant factors. This correspondence allows classical results from combinatorics and number theory to be applied directly to Boolean network dynamics. We find three distinct asymptotic scales for both average and maximum attractor lengths: typical networks scale as $\exp[(1/8+o(1))(\ln N)^2]$, the ensemble means grow as $\exp[N^{1/3+o(1)}]$, and extremal networks attain $\exp[(1+o(1))\sqrt{N\ln N}]$. Thus, ensemble averages are governed by rare network realizations and are unrepresentative of typical dynamics.

Biological systems are inherently multilayered: the same entities—genes, cells, or bacterial species—participate simultaneously in qualitatively distinct types of interactions, each carrying complementary information that no single relational view can capture. Analysing such systems with single-layer tools, or by collapsing layers into a monoplex projection, systematically discards inter-layer dependencies and can yield misleading conclusions about centrality, community structure, and robustness. The multilayer network formalism addresses, and muxViz established one of the first comprehensive toolkits for its structural analysis, but its R-only interface and dense data structures limit applicability to large biological networks. We introduce muxvizpy, a Python library that reimplements and extends the muxViz analytical catalogue with a sparse linear-algebra stack built on SciPy and PyTorch. Muxvizpy exposes seven categories through a unified, composable API and is numerically validated against muxViz on synthetic Erdős–Rényi and Barabási–Albert multiplex networks while substantially reducing peak memory and wall-clock time at scale. We illustrate its applicability on a virus–human protein-interaction multiplex in which computing some structural analysis was unfeasible. \\[2pt] muxvizpy is freely available under the MIT licence at https://github.com/CoMuNeLab/MuxVizPy. Mathematical definitions of all implemented metrics are provided in the Additional File.

In cell biology, computational models of biological systems range from minimal representations with a handful of parameters to whole-cell simulations tracking thousands of molecular species across a complete cell cycle. While these models span a continuum of detail, increasing complexity changes what they capture and are able to explain, what they can predict, and how they can fall short. A model contributes to understanding only when it makes novel predictions, reveals an unexpected coupling between processes, or fails in a way that identifies missing parameters. We contend that what is important for understanding is not the number of components or spatial dimensions a model contains, but the ratio of free parameters to the experimental constraints available to pin them down, and whether we can see why it produces the behaviors it does. Large-scale agent-based models of cytoskeletal dynamics or tissue mechanics that are built on a small number of physically grounded rules can reveal rich self-organization behavior precisely because their parameter spaces are small enough to explore systematically. By contrast, when free parameters grow faster than the data available to constrain them, models become progressively harder to interpret – and even disprove –regardless of their biological scope. We argue that the field needs to reconsider the goal of complex models. We should move away from trying to include as many parameters as possible and instead aim for systematic comparisons with simpler representations, dynamical analysis, and explicit model hierarchies that trace how cellular behavior emerges from its parts.

A central task in virtual cell modeling is predicting single-cell transcriptional responses to unseen genetic perturbations and drug combinations, and biological networks provide valuable priors on gene relationships. Existing graph-based models commonly use the same network to structure gene representations and mediate intergene interactions, thereby implicitly treating stable associations as perturbation-response pathways. Gene Ontology and control-derived coexpression networks encode relatively stable relationships rather than intervention-specific response directions or magnitudes. We therefore propose GeneGeoFlow, which conditions a control-anchored residual flow on gene-wise geometry derived from biological networks to learn intervention-specific transcriptional responses. GeneGeoFlow derives multi-scale spectral coordinates from Gene Ontology and control-derived coexpression networks. A perturbation-conditioned, gene-wise gating module selects relevant structural scales and network sources, yielding intervention-specific gene geometry. The resulting geometry conditions a control-anchored residual flow without explicitly propagating target-derived signals along the graph. Condition-wise optimal transport couples unpaired control and perturbed populations for training, while a Delta-correlation objective aligns the predicted and observed condition-level expression-shift directions. GeneGeoFlow achieves Pearson Delta scores of 0.8979 on the Norman additive benchmark and 0.9088 on five held-out drug combinations in the fixed ComboSciPlex test split. These results support perturbation-conditioned gene geometry as an effective structural prior for intervention-specific response prediction, without conflating stable gene relationships with response propagation.

Proteins occupy heterogeneous free-energy landscapes in which high-entropy ensembles converge toward compact, low-energy basins with multiple sub-states. Molecular dynamics can access these landscapes at atomic resolution, but exhaustive sampling remains computationally demanding. Meanwhile, most quantum approaches target only single optimal structures, leaving full ensemble energetic heterogeneity unexplored. We introduce a residue-level, gate-based quantum circuit framework for coarse-graining protein thermodynamics. Each amino acid is represented as a two-state qubit (stabilised vs. excited solvation state) based on residue solvation energetics. A structure-informed entanglement block then encodes covalent and non-covalent contacts using parameterised controlled gates, embedding correlations across the residue-interaction network. Sampling the circuit ($\sim 10^6$ measurements) yields binary thermodynamic microstates used to compute protein energy distributions, residue-level statistical couplings, energetic sensitivities, and information gains relative to total free energy. We showcase the framework on the benchmark Trp-cage miniprotein 1L2Y (TC5b) and 9GDL, a disulfide-stabilised Trp-cage-fortified exenatide chimera. For 1L2Y, the circuit reproduces a structured, folding-funnel-like energy distribution. Comparative analysis with 9GDL reveals shifts in global energy distributions and residue-level stability profiles. Coupling and information-theoretic analyses localise residues associated with ensemble reorganisation, while multi-body couplings show the circuit resolves both direct and indirect statistical correlations. This framework expands quantum protein modelling beyond single-structure optimisation toward ensemble-level characterisation, capturing key features of rugged energy landscapes to guide protein design, mutation mapping, and allosteric pathway identification.

Excitable and oscillatory responses are widespread in cell signaling, but how these arise in the same system is unclear. We propose a minimal model of a saddle node bifurcation on an invariant circle, where slow fluctuations modulate transitions between excitable and oscillatory states, and noise drives excitability. Fitting the model to cell signaling data from embryonic stem cells reveals that these cells operate close to the bifurcation. Proximity to the bifurcation allows the cell to tune different dynamic states, suggesting a general principle for versatile cell signaling.

Abstract reaction networks appear not only as models of chemical reactions but also as models of complex systems, with applications in areas such as ecology and epidemiology, and as one of several alternative paradigms for non-standard computation. It is therefore of interest to determine whether an abstract reaction network can be realized by a concrete set of molecules and plausible chemical reaction mechanisms. It is known that a reaction network has a realization in terms of chemical graphs (i.e., Lewis structures) if and only if it is conservative. Here we consider the problem of assigning a set M of known molecules to a set X of abstract entities in a reaction network (X,R) such that each reaction satisfies mass balance and adheres to one of an allowed set of chemical reaction mechanisms. We show that this problem is NP-complete. Nevertheless, it can be solved in practice using a backtrack-and-prune algorithm inspired by the VF2 family of algorithms originally designed for the subgraph isomorphism problem.

Fine-tuning a large language model on new data degrades what it previously learned. We present Omega-S, a drop-in penalty computed from the weight matrix alone: it needs no previous-task data, no Fisher matrix and no stored copy of the old weights. It is three lines in an existing training loop and adds under 4% to the cost of a step. Retention. On Llama-3-8B with LoRA, fine-tuned from code to prose and measured by HumanEval over ten seeds, Omega-S retains more of the original capability than no regularisation on 9 of 10 seeds (0.173 -> 0.238 absolute pass@1; sign test one-sided p=0.011, Wilcoxon p=0.006), as a retention ratio, 62.9% -> 84.1%. It also beats tuned weight decay on 10 of 10 seeds (p=0.002) and tuned EWC on 8 of 10 (p=0.014), every arm re-measured in the same session. Mechanism, measured rather than asserted. Omega-S is topological by construction, its objective built from Tr(A^3), but we measured which of its four factors actually moves and three do not: their elasticity with respect to the weights is at or below 1e-4, against 9e-3 for the degree-variance term. As implemented, the composite reduces to a penalty on the variance of node degrees, which means row magnitude in square modules and directional alignment in non-square ones. We report this because a method whose name promises one thing and whose gradient does another should say so. We also enumerate the open design choices, including a contrast-preserving construction that does what it was designed to do and makes retention worse on all ten seeds. Repeating an identical configuration, same seed and same hardware, gives a standard deviation of 0.104 in retention ratio. We have not found this quantified for low-rank fine-tuning of language models, and it bounds every seed-paired comparison in this literature, ours included. Code, per-seed results and the full record of negative results are available.

Identifying structural features ensuring the persistence of a reaction network is fundamental for understanding the evolution and complexification of biological entities. By defining organizations as subsets of species that are closed and self-maintaining, Chemical Organization Theory (COT) shows that organizations enable filtering out the regions of the phase space where attractors can exist. Despite various theoretical explorations current methods to compute the organizations face unclear combinatorial challenges. In this article, we first perform a systematic study of such combinatorial challenges and identify a criteria of productive novelty and irreducibility that separate relevant from redundant combinations. Second, we identify the minimal building blocks that shall be combined to build organizations, called elementary reaction closures (ERC), and characterize them as a hierarchy. Third, we show that all persistent modules can be generated as combination of ERCs only, and operationalize two properties among ERCs, synergy and complementarity, that define a sufficient criteria to build such generators in a minimal way. Fourth we show that every relevant persistent module can be built from such minimal sequences of ERCs. We next show that the synergies and complementarities we use to generate persistent modules not only scale radically slower than usual combinatorial methods. We quantify these structures across $438$ biological reaction networks from the BioModels and BiGG databases. The proportion of fundamental synergies and complementarities among all ERC pairs shrinks as networks grow even as raw counts increase, confirming that larger biological networks achieve persistence through a progressively smaller and more selective sets. We do not prescribe algorithms to compute organizations, discuss how algorithms based on our results could permit applying COT.