Papers

A central hypothesis in aging biology proposes that the gradual accumulation of stochastic molecular damage drives transcriptional deregulation (TD), ultimately leading to functional decline. However, existing single-cell analyses largely rely on static expression measurements and provide limited resolution of how regulatory responses change at the level of individual cells. Here, we leverage DKO, a deep learning framework that predicts cell-specific virtual gene perturbation effects from unperturbed transcriptomes, capturing cell-to-cell differences in regulatory context without requiring a reconstructed gene regulatory network or matched perturbation data. Applying DKO to single-cell transcriptomic datasets from humans, mice, Drosophila, and C. elegans revealed that the transcriptional consequences of perturbing the same gene became increasingly heterogeneous across cells with age, indicating that gene regulatory responses become less consistent and more dependent on cellular context. Furthermore, we find that cells with higher somatic mutation burdens have greater knockout-response variability, linking accumulated molecular damage to regulatory inconsistency. Within the same cellular populations, higher impact variability was also associated with higher senescence scores, linking transcriptional dysregulation to senescence-associated states. Our findings identify a progressive loss of regulatory robustness at both the gene and cell levels as a conserved feature of aging, linking stochastic molecular damage to increasing heterogeneity in cellular regulatory responses.

We investigate the influence of decay rates on the dynamics of three-dimensional gene regulatory networks modeled by nonlinear differential equations with logistic sigmoidal regulatory functions. Particular attention is paid to the existence, location, and stability of equilibrium points in systems involving numerous interacting parameters. An explicit analytical representation of a central equilibrium is obtained under suitable relations between the regulatory weights, thresholds, and decay coefficients, allowing its location to be determined without solving the corresponding nonlinear algebraic system numerically. General formulas for the characteristic polynomial of the linearized system are derived, and spectral conditions for instability are established for several classes of regulatory interactions. We demonstrate how variations in decay rates can change the number and stability of equilibria, transform oscillatory dynamics, and influence the emergence or suppression of complicated attractors. The possibility of chaotic behavior is investigated through eigenvalue configurations, Lyapunov exponents, and the Kaplan-Yorke dimension. Numerical examples illustrate the sensitivity of the dynamics to decay parameters, including the suppression of previously observed chaotic behavior. The results provide an analytical framework for investigating equilibrium configurations, stability transitions, and complex dynamics in parametrized gene regulatory networks.

Subcellular signaling pathways regulating, e.g., cell division or synaptic plasticity, are often modeled as dynamical networks of reactions satisfying the law of mass action. To remain physically grounded, these models must satisfy different physical conservation laws, such as conservation of energy. This gives rise to thermodynamic constraints formalized by the Wegscheider conditions, which ensure that the system is in detailed balance at thermodynamic equilibrium. Explicitly imposing Wegscheider conditions requires identifying cycles in the reaction network. This can become computationally cumbersome, especially in large, complex systems. Furthermore, even when these constraints are explicitly available, enforcing them requires choosing a reparameterization of the network, which may introduce asymmetries among the parameters and affect numerical conditioning. Moreover, this approach often causes inverse uncertainty quantification methods like Markov chain Monte Carlo (MCMC) to struggle. Instead of imposing Wegscheider conditions directly in the parameter space, we propose an approach for promoting physical consistency by penalizing parameters that violate detailed balance. We first illustrate the approach on a simple reaction network containing a single cycle. We show that parameter values inferred without thermodynamic constraints can reproduce the experimental data seemingly accurately, while substantially violating the Wegscheider condition. In contrast, the proposed approach preserves the agreement with the observed dynamics while yielding posterior samples that approximately satisfy the thermodynamic constraint. We further apply the approach to a model of PKA-calcineurin signaling in synaptic long-term depression, demonstrating its applicability to a more complex biochemical reaction network.

Deuterium metabolic imaging (DMI) measures glucose-derived labeling in vivo, but noisy time courses complicate network-level interpretation. We developed an uncertainty-aware workflow that integrates dynamic DMI data with genome-scale metabolic modeling (GEM) through soft flux constraints. After injection of [6,6'-2H2]glucose, dynamic spectra were acquired from gastrocnemius, soleus, and tibialis anterior muscles of control and high-fat diet (HFD) mice. A hierarchical four-state kinetic model estimated labeling parameters while accounting for shared plasma input, tissue-specific kinetics, voxel-level fluxes, isotope effects, and measurement uncertainty. Because glucose labeling was more reliable than lactate and glutamate/glutamine signals, only two DMI-derived quantities were used as GEM anchors: glucose processing mapped to hexokinase flux, and glucose-derived oxidative entry mapped to mitochondrial pyruvate dehydrogenase flux. Two-sided soft target intervals with weighted penalties defined DMI-compatible feasible regions for each mouse and muscle. Flux variability analysis and sampling assessed changes across reactions and metabolic subsystems. The constraints were compatible with the muscle network and affected feasible activity in glycolysis, lactate metabolism, pyruvate oxidation, the TCA cycle, oxidative phosphorylation, and redox pathways. HFD-associated changes differed by muscle: tibialis anterior and soleus showed reduced glucose-linked feasible activity, whereas gastrocnemius showed mixed changes with preserved or increased TCA and redox activity. This framework connects noninvasive tracer imaging with network-level metabolic inference while avoiding strong claims from poorly identified signals.

This study aimed to reconstruct a regulatory network suitable for systems biology analysis of diabetic bone marrow. Three artificial intelligence large language models independently analysed 18 Discovery and 43 Validation publications. Investigator guided integration produced a frozen network of 62 biological entities and 79 directed interactions. Open source software Gephi created a visual representation before and after main node removal. The plausibility of the reconstructed interactions was tested against independent curated knowledge resources, including the Kyoto Encyclopedia of Genes and Genomes, connectomeDB2025, Reactome, and Gene Ontology. Further analysis examined physiological correspondence, missing edge sensitivity, and graph constrained interrogation. Functional modules encompassed endothelial perivascular support, stem cell retention and trafficking, neural circadian regulation, inflammatory myelopoiesis, metabolic stress, and epigenetic regulation. They intersected through shared effects on endothelial function and stem cell behaviour. The endothelium was the principal connecting node. Computational removal of the endothelial node increased network fragmentation and diameter, while average directed path length remained essentially unchanged. External resources supported 60/79 interactions (75.9%); 19 remained unresolved, with none classified as contradictory. Complementary assessments supported biological plausibility and faithful interrogation of the reconstructed network. This investigator led dialogue with artificial intelligence provides a foundation for a systems level understanding of diabetic bone marrow dysfunction. The approach may facilitate the generation of testable hypotheses and the identification of mechanistic targets for coordinated restoration of disrupted niche interactions.

We introduce positively transverse reaction networks and show that, with mass-action kinetics, no bounded trajectory with positive initial condition has a boundary omega-limit point. These networks are characterized by a geometric property of a collection of polytopes defined by a subset of complexes of the network.

Repairing scientific models from functional observations differs fundamentally from supervised prediction: feedback may certify a solution without revealing which structural correction is responsible. We study this setting for genome-scale metabolic model (GEM) repair, where multiple reaction edits can explain the same phenotypes and many apparently distinct edits correspond to the same biological mechanism. This many-to-one structure creates a hidden failure mode for conventional exploration: diversity in the output space need not translate into diversity of scientific hypotheses. We introduce QuotientPO, which collapses equivalent repairs into canonical mechanisms and optimizes exploration directly over the resulting quotient space. To make quotient exploration informative under finite rollouts, we derive a kernelized Rényi estimator that resolves graded crowding among distinct repair cores beyond coarse exact-match counts. On 2,212 held-out GEMs, QuotientPO improves Success@32 from 17.93% to 20.10% (+12.1% relative) while consistently increasing distinct successful-core discovery under the same sampling budget. These results establish quotient-space exploration as a principled approach to mechanism-level discovery under verifier-induced equivalence.

Ayurvedic formulations represent a rich source of traditional knowledge, encompassing distinct therapeutic approaches, medicinal ingredients, and underlying principles of drug action. However, systematic integration of this knowledge with standardized modern biomedical terminology and resources remains limited. Here, we present IMPPAT 3.0, an updated version of the Indian Medicinal Plants, Phytochemistry And Therapeutics (IMPPAT: https://cb.imsc.res.in/imppat/) database that extends its scope by incorporating curated and harmonized information on Ayurvedic formulations from authoritative sources, while expanding its coverage of medicinal plants, phytochemicals, and therapeutic use terms. IMPPAT 3.0 further enriches plant and phytochemical information with annotations on clinical trials, patents, bioactivity, protein targets, spectral data, plant image links, and AyuCaRe reports. The compiled data is also organized as IMPPAT-KG, a knowledge graph that enables exploration of relationships among IMPPAT entities. We demonstrate the utility of the resource through a network medicine analysis of single herbal drugs (SHDs) reported for rheumatoid arthritis (RA), revealing significant proximity of their protein targets to RA-associated genes, and pinpointing of phytochemical pairs with potentially complementary effects on RA pathogenesis. In conclusion, by integrating traditional medicinal knowledge, IMPPAT provides a FAIR-compliant ethnopharmacological platform for systematic investigation of Indian medicinal plants and Ayurvedic formulations for drug discovery applications.

The dynamics of chemical reaction systems with two chemical species can be described (under mass-action kinetics) by planar autonomous systems of ordinary differential equations (ODEs) with right-hand sides containing polynomials. While their behaviour close to hyperbolic critical points is well understood, their dynamics has not been fully characterized for strongly non-hyperbolic critical points, where the linearization close to the equilibrium is given by the zero matrix. In this paper, the behaviour of such chemical reaction systems is investigated. Considering second-order chemical reaction networks with a positive equilibrium point, we show that only nine geometric equivalence classes are possible out of the sixteen classes for general ODE systems. Examples of chemical reaction networks in each class are presented. All sixteen geometric equivalence classes can be realized by chemical systems if the non-hyperbolic equilibrium point is located at the origin.

Identifying kinetic parameters and conserved quantities that support multistationarity and multistability remains challenging in chemical reaction network analysis. Structural theory provides existence and exclusion results, but does not always yield explicit, interpretable descriptions of the corresponding parameter regimes. We present a computational and machine-learning framework for discovering and characterizing parameter regimes with numerical evidence of multistationarity and multistability in mass-action systems. We refer to these numerically detected behaviors as approximate multistationarity and approximate multistability. The framework combines Latin hypercube sampling, adaptive sampling, and interpretable surrogate modeling. Logistic-regression and neural-network surrogates are trained separately on accumulated detection labels and independently guide later sampling toward regions with high scores for approximate multistationarity or multistability. Boundary refinement targets transitions between detected and undetected numerical classes. Buckingham–Pi groups are then used with CN2 rule induction to obtain interpretable empirical descriptions of detected regimes. Illustrative networks yield rules consistent with known multistationarity conditions, while the Schlögl model demonstrates adaptive enrichment of multistationarity detections. In the hybrid histidine kinase model and the sequential dual phosphorylation cycle, adaptive sampling substantially increased the number of parameter - total pairs with detected multistability relative to global sampling under equal sampling budgets. The framework complements analytical chemical reaction network theory by providing efficient, parameter-resolved discovery and characterization of multistationary and multistable regimes.

Protein function annotation needs to know which predictions to distrust, not only what a model predicts. We ask whether tissue-specific interaction structure carries that information. Our candidate signal is effective resistance, used previously to relieve over-squashing by rewiring. Across 24 tissue-specific interactomes it is dominated by inverse degree, and the degeneration deepens as the co-expression filtered network grows, with a Spearman correlation of -0.955. The residual departure from that limit exceeds degree-preserving null graphs in all 24 networks. Controlling for predictive entropy, degree, annotation cardinality, local structure and feature-only difficulty, the residual explains additional per-node loss in 19 of 24 held-out networks once a permutation floor is subtracted, at every depth, and the effect strengthens monotonically with depth. The increment reaches 0.37% of the variance the controls leave unexplained, 5.6 times a permutation floor, against 1.5 times when the model is retrained in a degree-preserving null world. Selective prediction improves negligibly. The signal is reproducible; degree degeneration bounds it.

The next-generation matrix method is a powerful tool for computing the basic reproduction number in compartmental mathematical models of infectious diseases. The method has been recently extended to mathematical biochemistry, where it can be used to establish parameter regions of stability and instability for boundary steady states. Several significant challenges in its application remain, however, particularly around establishing conditions under which the method is mathematically valid, computationally tractable, and biologically meaningful in the biochemical setting. In this paper, we address these challenges by shifting the interpretation from new infections in the epidemiological setting to autocatalysis in the biochemical setting. We introduce a graph, called the ALT-graph (autocatalysis-leak-transition graph), which decomposes the contribution of each reaction to the Jacobian as autocatalytic, leak, or transition edges. We then present a systematic algorithm for splitting the ALT-graph, which is guaranteed to produce a valid reproduction number, $\rho(FV^{-1})$, while also decreasing computational complexity by lowering the rank of $FV^{-1}$. We apply the method to models of both biochemical reaction networks and infectious disease spread.

Phenotypic outcomes arise from coordinated changes across molecular layers. Although multi-omics data capture complementary molecular signals, integrating them into coherent biological mechanisms remains challenging due to limited overlap at the gene level. Network propagation maps significant gene hits from each dataset onto a prior knowledge molecular interaction network and diffuses their signals to neighboring genes. Mathematical integration of propagated networks, followed by clustering and pathway enrichment, reveals shared molecular networks, pinpointing the core biology captured across orthogonal molecular layers.

Biological models based on mathematical simulations (e.g., ordinary differential equations) are hugely beneficial for understanding biological processes and predicting the outcomes of experiments. Such models are built using a variety of sophisticated modeling tools. Unfortunately, the complexity of these tools is often a barrier to use by experimental biologists. To improve the usability and impact of biological models, we have created a desktop platform (DeskIridium) and a web-based simulator (WebIridium, hosted as a GitHub page) that, in combination, provide features that promote ease of use by experimental biologists, including: human-readable editing of models via the Antimony syntax; sliders for adjusting parameter values while simultaneously displaying simulation results in real time; zero-install distribution; AI chat integration; and generation of standalone Python implementations of SBML models to support reproducibility. The desktop version uses the well-established libroadrunner simulation package; the web-based application uses an Emscripten translated version of COPASI. Both simulators are SBML-compatible, allowing users to easily import previously built models or edit existing models using the easy-to-use antimony language. We used modern AI methods for software development and share the lessons we learned. Binaries, source code, and the web interface are available at: https://github.com/sys-bio/IridiumSimulator and https://github.com/sys-bio/WebIridium.

Genetic design automation maps Boolean specifications to regulatory networks and DNA sequences, but a successful mapping does not establish the probability that every output remains correct over time. This paper presents Proof-Carrying Synthetic Biology (PCS-Bio), a formal architecture for acyclic genetic logic circuits operated under fixed inputs. An ideal checker connects Boolean equivalence, sequence and model provenance, local stochastic contracts, and any surrogate-to-target discrepancy bounds. Under explicit semantic and checker assumptions, these obligations imply a joint temporal refinement guarantee without assuming independence among molecular events. Assignment-specific sufficient cubes select the branches needed to prove each output value. This can reduce risk and delay accounting provided that the retained local contracts remain uniform over the physically present but logically masked inputs. The reported evaluation exercises three separate interfaces: an archived Cello execution with sequence reconstruction, a restricted exact rational checker for a synthetic finite-state instance, and diagnostic IEEE~754 binary64 replay. Across 60 deterministic Boolean formula graphs, slicing reduced the worst-assignment risk-charge sum in 55 cases and the maximum accounting delay in 54. These results concern proof accounting and selected implementation components. They do not establish an integrated sequence-bound temporal certificate or physical validation.

Growing biochemical systems are intrinsically noisy, with one important source of stochasticity arising from the discrete reaction events of the underlying biochemical network. Because growing systems do not generally admit stationary abundance distributions, it is unclear how to separate this intrinsic chemical noise from other sources of variability in population and single-cell data. Here we develop a large-deviation theory for exponentially growing stochastic chemical reaction networks by decomposing abundance into volume and composition. In these coordinates, balanced growth corresponds to a stable composition together with exponentially increasing volume. We show that composition fluctuations satisfy a large-deviation principle with speed equal to the volume of the growing system, and the corresponding quasipotential is selected by a solution of the contact Hamilton–Jacobi equation. We also derive a fluctuation theory for accumulated observables and show that their covariances decay with accumulated-volume, rather than with physical time. Applications to a minimal autocatalytic network and a coarse-grained cellular growth model demonstrate how the theory predicts composition quasipotentials, growth-rate fluctuations, and correlations between observables. The framework therefore provides a stochastic null model for single-cell heterogeneity generated by the intrinsic stochasticity of the underlying chemical reaction network.

We prove the Global Attractor Conjecture for complex balanced mass-action reaction networks whose reachable siphons satisfy two structural conditions, allowing multiple linkage classes. The proof proceeds in two steps. A structural condition relating the stoichiometric space to the reactions active within a boundary face ensures that a stoichiometric compatibility class contains at most one boundary equilibrium with any prescribed zero set. Finiteness of the possible zero sets and connectedness of the $\omega$-limit set then imply that any boundary limit set consists of a single equilibrium. To exclude convergence to such an equilibrium, we consider the embedded reaction network obtained by projecting onto the vanishing species and freezing the concentrations of the surviving species at a positive limit. The structural condition guarantees complex balance of this embedded network, while a further condition on its minimal active linkage classes ensures that an explicit Chetaev function is strictly increasing near the boundary. This excludes the boundary point as an accumulation point whenever the surviving concentrations converge. Consequently, every positive trajectory converges to the unique positive equilibrium in its stoichiometric compatibility class. Examples illustrate the hypotheses and their relation to strong endotacticity.

Autocatalysis, the ability of a chemical subsystem to sustain its own constituents when supplied with sufficient food molecules, has been closely related to the origin of life on Earth. Emerging from Wilhelm Ostwald's considerations about an explicit autocatalytic reaction, different notions of autocatalysis have been developed over the years. The two most prominent are reflexively autocatalytic F-generated sets (RAFs) and stoichiometric autocatalysis. After having shown that each RAFs is, under reasonable conditions, in general stoichiometrically autocatalytic, we examine here the relationship between the two notions of minimality: irreducible RAFs and autocatalytic cores. To this end, we overcome the obstacle that RAFs and stoichiometric autocatalysis have been formalized in distinct systems of chemical reactions, i.e., catalytic reaction systems (CRS) and chemical reaction networks (CRNs), respectively. We show that reactions in a CRS constitute equivalence classes of reactions of the corresponding CRNs w.r.t. their specific catalyzations. Using the fact that CRN and CRS can be canonically identified whenever each CRS reaction is associated with a single catalyzation, we demonstrate that the Kőnig graph of a monocatalyzed, irreducible RAF is composed of strong blocks devoid of food and waste species that are separated by reaction vertices, each of which contains an autocatalytic core. In fact, a single irreducible RAF can, in general, contain exponentially many autocatalytic cores.

Living systems compute with chemistry by mapping environmental signals onto specific internal chemical states. Despite recent advances in molecular programming, it remains unclear which physicochemical features control the computational expressivity of chemical systems. Here we inverse-design thermodynamically consistent chemical reaction networks whose steady-state response to an environmental input computes a target nonlinear function. Using implicit differentiation we train the free-energy landscape directly: standard chemical potentials, transition-state energies and thermodynamic drives. Increasingly large networks generated by elementary ligation and cleavage steps fit increasingly complex nonmonotonic polynomial functions, with expressivity scaling logarithmically with network size, predicted primarily by the number of reactions. Training individual energetic parameter classes reveals that internal thermodynamic drives, capable of breaking detailed balance, dominate trainability, with comparable performances achieved only by pairs of parameter classes. These results identify nonequilibrium drive as the most effective single resource for steady-state computational expressivity in chemical reaction networks.

Protein-protein interaction databases record evidence without weighing it. In BioGRID, an interaction asserted once by a single high-throughput screen and one confirmed by twenty laboratories across a dozen assays are the same kind of row in the same file. Tools built on such databases inherit that flattening: they draw every reported interaction as an edge, and the resulting picture states that two proteins interact without stating how much anyone should believe it. We present ProLiVis 2.0, a rewrite of the literature-centric visualization system of arXiv:2111.12794. It contributes three things. First, a citation-trust model that scores each interaction from seven terms, including a term for the number of independent laboratories behind the supporting publications, obtained by clustering those publications over shared institutional affiliations; a plain count of publications cannot distinguish five confirmations from one group publishing five times. Second, a deterministic reformulation of the center layout, closed-form and $O(n \log n)$, which replaces the force-directed placement of the original and makes published figures regenerable from a session manifest. Third, an implementation that runs entirely in a web browser, with an embedded analytical database, requiring no installation and uploading no data. On BioGRID release 5.0.260 restricted to SARS-CoV-2, 24,344 of 34,540 reported interactions (70%) rest on a single publication, and raising the trust threshold to 0.2 leaves 11,320 of them. That the large majority of a curated interaction network is unreplicated is a fact no existing view of the database makes visible.

The Interconnectedness Coefficient (IC) is a bounded semi-local graph-theoretic node measure designed to identify connector vertices between cohesive network regions. Such connector vertices, also referred to as bridging nodes, may mediate between locally cohesive regions even when they are neither hubs nor themselves highly clustered. The IC preferentially assigns high values to weakly clustered focal vertices whose adjacent vertices remain strongly clustered after exclusion of the focal connection. Candidate vertices are required to have degree at least two. The construction is partition-free, uses information within radius two, and requires no predefined community or module partition. The range and extremal properties of the score are derived analytically. Exact graph families isolate its maximal response to fully cohesive branches, its controlled response to a single cohesion defect, its invariance under a cohesion-free hub extension, and a sharp fragmentation threshold. A separate application to a Human Interactome Map reveals pronounced degree-dependent stabilization of IC values near the network's mean clustering level. This behavior follows directly from the multiplicative definition. If focal clustering tends to zero while mean leave-one-out cohesion in the neighborhood stabilizes, the IC converges to that neighborhood-cohesion level. Among the highly ranked IC vertices are proteins with established interface, scaffold, and adaptor roles in molecular complexes. The IC is therefore positioned as a semi-local connector measure for cohesive network regions.

Protein-protein interaction (PPI) databases do not faithfully reflect biological realities. Instead, they are influenced by study and technical biases that distort certain protein and interaction attributes. Machine learning models can exploit these as learning shortcuts if the negative dataset is not constructed with care. So far, the shortcuts introduced during PPI dataset construction have only been examined in isolation. Here, we systematically characterize both reported and, to our knowledge, previously unreported biases in PPI datasets that lead machine learning models to learn shortcuts instead of biological signal. We analyze HIPPIE, IntAct, and STRING, dedicated PPI databases, as well as two datasets derived from 3D-structural information in the Protein Data Bank (PDB). We show that random data splitting introduces strong topological shortcuts. When train-test protein overlap is removed, the resulting datasets still retain usable shortcuts stemming from self-interactions, taxonomic identity, and functional relatedness, whose prevalence interestingly depends on the data source. We further show that sampling negatives from a set of high-confidence non-interactors, an intuitively appealing choice, can amplify the shortcut stemming from functional relatedness. To detect and mitigate these biases, we provide an open Nextflow pipeline that combines similarity-aware, data-loss-minimizing dataset splitting with bias-minimizing negative sampling, both formulated as integer linear programs. Its key concept of quantifying biases to minimize them through optimization-based negative sampling can, in principle, be extended to any machine learning problem where the pool of negative candidates is much larger than the positives and is thus of interest also beyond PPI prediction.

Stochastic switching between gene expression states, coupled with production and degradation dynamics, governs the accumulation of mRNA and proteins in cells. The concentrations of these accumulated entities dictate the phenotypic distribution of genetically identical cells. The underlying accumulation dynamics are well-captured by a two-state promoter switching model, with statistical and thermodynamic properties quantified via the Fano factor and entropy production rates. However, how these measures correlate with concentration distributions and their shifts under varying kinetic parameters remains largely unexplored. To this end, we use chemical master equations to study a generalized model of mRNA accumulation dynamics in the presence of stochastic switching between two activity states and state-dependent production and degradation rates. We derive exact expressions for the steady-state probability distribution and analytically compute the mean concentration, Fano factor, and entropy production rate (EPR). Simplifying these expressions, we identify contributions arising from stochastic switching rates and relaxation dynamics toward equilibrium in each activity state. Next, using our theoretical results, we characterize the variation in the Fano factor and EPR as a function of mean expression during modality changes of the distributions mediated by the variation of switching rates. We also identify the conditions in kinetic parameters that achieve the highest Fano factor and entropy production rates. Our findings establish a generalized framework for examining stochastic accumulation dynamics, clarifying how kinetic parameters dictate molecular distributions, noise, and dissipation. These insights extend readily to broader contexts coupling stochastic switching with accumulation, including protein burst dynamics, phenotype-switching-mediated drug intake, and queuing theory.

Protein-protein interaction (PPI) network alignment combines topological and sequence information to identify conserved modules across species, but global alignment remains challenging: heuristics sacrifice optimality, while exact methods lack scalability. We model the alignment as a weighted maximum common induced subgraph problem and reformulate it through the modular product graph to a minimum-weight vertex cover on the complement, with node weights carrying sequence similarity. To solve this problem, we develop a hybrid framework combining kernelisation, branch-and-bound, and seven Quantum Approximate Optimisation Algorithm (QAOA) formulations. These formulations differ in how the cover constraints are enforced, from penalty terms in the cost Hamiltonian to mixers confined to the feasible subspace. For single round QAOA, we derive closed-form expressions for the expected cost of four circulant mixer variants, enabling performance characterisation without circuit simulation. Applied to synthetic and real-world networks reduced to KEGG pathways, the QAOA formulations achieve high topological conservation on the aligned core while at least maintaining biological conservation comparable to leading classical aligners, at the cost of reduced node coverage. Across selected KEGG pathways, the aligned subnetworks retain disease-associated proteins, preserving biologically relevant information. Cheaper formulations leave more edges uncovered, while enforcing feasibility in the mixer raises circuit depth by one to two orders of magnitude. Together, these results highlight the potential of quantum optimisation for PPI network alignment and the resource trade-offs that will shape its scalability as quantum hardware matures.

The purpose of this work is twofold. In the first part, we consider layered mixed matrices introduced by Murota, relate them to existing notions in combinatorial commutative algebra, and investigate the irreducibility of their determinants. Furthermore, for a layered mixed matrix in combinatorial canonical form, we determine the sparsity structure of its inverse. That is, we characterize which entries of the inverse are nonzero. In the second part, we establish for the first time a formal connection between these algebraic results and the theory of buffering structures for reaction networks developed by Mochizuki and Okada. We identify the lattice of buffering structures with the lattice of order ideals of the block poset of the combinatorial canonical form of the associated layered mixed matrix. This allows us to characterize the reducibility of the symbolic Jacobian determinant as a polynomial in the reaction-rate derivatives, as well as the nonzero sensitivity responses of species concentrations to reaction-rate perturbations.

Many natural and engineered systems are large complex networks of interacting components, and external perturbations drive them along different dynamical paths. Identifying which perturbation occurred matters for diagnosis, control, and prediction. Yet often times only a few components can be jointly monitored. Which components should be monitored? Experimental practice usually favors placing reporters at the most sensitive sites, where perturbations produce the largest effects. Using a simple dynamical model for complex systems with heterogeneous connectivity, we ask how sparse reporter panels should be chosen to classify shocks from partial trajectories when repeated trials only approximately reproduce an ideal initial condition. Once that reproduction is imperfect, sensitivity ranked panels fall far short of optimal, and the shortfall grows with the noise. We find that the best panels mix two kinds of reporters. A promiscuous reporter responds to most shocks, so it separates them mainly by degree, and degree fluctuates from trial to trial. A dormant reporter responds to only a few shocks but strongly, and its answers do not scatter as much between trials. As noise grows, the cost of losing a dormant reporter rises to meet the cost of losing a promiscuous one. Panels of either kind alone classify worse than the mixture, and no property of the members collected individually explains the ordering. Most of all, we find that only a minuscule number of reporters are needed on a panel to accurately identify which shock hit the system. We implement an efficient algorithm to assemble such a panel. Together these results provide a low cost practical design principle for monitoring large complex dynamical systems.

Protein co-abundance clusters preserved across tissues can reveal shared disease mechanisms and candidate therapeutic targets, particularly when proteins implicated in organ-confined diseases converge in peripheral or accessible tissues. However, previous cross-tissue studies have focused on biologically pre-selected tissue pairs, leaving most possible combinations and non-obvious relationships unexplored. We present an LLM-agent framework for large-scale, evidence-grounded comparison of tissue-specific protein co-abundance networks. The framework constructs tissue networks, derives pairwise consensus clusters, and integrates evidence from expression atlases, protein interaction and complex databases, pathway annotations, disease catalogues, and literature. Applied to all 820 pairwise combinations of 41 human tissues and fluids, it identified 1,833 conserved co-abundance clusters across 406 tissue pairs. Colon, synovial fluid, blood, cerebrospinal fluid, and bone marrow were the most broadly connected tissues, while the most cluster-rich pairs were dominated by bone marrow. The analysis also highlighted non-obvious relationships: skin-bone marrow exceeded the anatomically adjacent bone-bone marrow pair, while colon-breast contained cancer-relevant clusters involving extracellular-matrix remodeling, lipid metabolism, and immune modulation. Cluster-level analyses generated further mechanistic hypotheses, including a brain-gut extracellular-vesicle/redox/serotonin-cofactor axis and a liver-bone marrow stress-response axis involving genes linked to white matter disease. These results provide a global, comparable landscape of conserved protein co-abundance and a hypothesis-generating resource for mechanistic and therapeutic exploration. Code and data are available at https://github.com/Gry1005/AgenticAI-conserved-cross-tissue-protein-co-abundance.

Orchestra composes two independently built bioinformatics MCP servers – RegNetAgents, which infers gene regulatory network topology from ARACNe networks, and CASCADE, which supplies four independent evidence sources (LINCS knockdown, DepMap essentiality, super-enhancer status, DoRothEA transcription-factor confidence) – into one multi-agent workflow exposed via the Model Context Protocol. Its central architectural claim is that requiring RegNetAgents' topology evidence and CASCADE's experimental evidence to agree on a candidate regulator yields a more trustworthy candidate than either alone – not previously tested directly, since RegNetAgents' own validation asked only whether its candidate lists beat chance. We test this on the TCGA tumor-acquired regulator tier (regulators in a gene's tumor ARACNe network but absent from the GREmLN population-averaged baseline), selecting candidates by ARACNe mutual-information (MI) edge weight. On RegNetAgents' published BRCA/COAD focal-gene panel plus matched negative controls, agreement among at least 2 of the 4 CASCADE sources predicts OncoKB cancer-gene status among focal genes (odds ratio 2.89, Benjamini-Hochberg-adjusted p=0.0166) but not among negative controls (p=0.0721); a single source is not diagnostic for either group. The pattern replicates and strengthens in a third cancer type, STAD, on a separately constructed panel (odds ratio 5.82), and against an independently curated ground truth (the Sanger COSMIC Cancer Gene Census). MI edge weight is the strongest single predictor overall (p=0.0003); a logistic-regression likelihood-ratio test confirms corroboration adds value beyond it in both panels (p=0.0234; p=0.0001). Every experiment invokes Orchestra's real agentic entry point.

Biological reaction networks often exhibit complex transient behavior that cannot be explained solely by the analysis of steady states or long-term persistence. Existing structural approaches identify persistent system properties and have considered transitions between organizations, but do not provide a general criterion for admissible transitions between arbitrary species configurations. We introduce admissible transitions, a new structural concept that describes feasible changes between species subsets using only reaction network structure and feasible reaction fluxes, independently of kinetic parameters. We prove that solutions of reaction-based ordinary differential equation systems induce canonical sequences of admissible transitions with a fundamental asymmetry: closure-building transitions are uniquely determined by network structure, whereas downward transitions are generally non-unique and depend on the realized trajectory. This establishes a structural layer linking network topology to transient system evolution. The framework is illustrated using a classical HIV immune-response model. By extending structural reaction network analysis from persistent states to transient dynamics, the proposed theory provides a general, kinetics-independent framework for analyzing reachability, organization formation, and transient behavior in biological systems.

This paper establishes that the system of ordinary differential equations arising from the sequential and distributive dual futile cycle has the structural capacity for Hopf bifurcations, whenever it is endowed with general parameter-rich kinetics, but it loses such capacity if it is endowed with mass action kinetics. The proof of the latter fact relies on a Routh-Hurwitz approach, improved by few preliminary structural considerations, but a decisive contribution came from ChatGPT Sol 5.6, which provided a nontrivial positivity certificate. A second purpose of the paper is indeed to document such AI-contribution, exploiting this well-studied toy-model, and to offer an example of the utilizability of such tools in the context of computer algebra and positivity certificates for large polynomials.