Search

While noise is generally associated with uncertainties and often has a negative connotation in engineering, living organisms have evolved to adapt to (and even exploit) such uncertainty to ensure the survival of a species or implement certain functions that would have been difficult or even impossible otherwise. In this article, we review the role and impact of noise in systems and synthetic biology, with a particular emphasis on its role in the genetic control of biological systems, an area we refer to as Cybergenetics. The main modeling paradigm is that of stochastic reaction networks, whose applicability goes beyond biology, as these networks can represent any population dynamics system, including ecological, epidemiological, and opinion dynamics networks. We review different ways to mathematically represent these systems, and we notably argue that the concept of ergodicity presents a particularly suitable way to characterize their stability. We then discuss noise-induced properties and show that noise can be both an asset and a nuisance in this setting. Finally, we discuss recent results on (stochastic) Cybergenetics and explore their relationships to noise. Along the way, we detail the different technical and biological constraints that need to be respected when designing synthetic biological circuits. Finally, we discuss the concepts, problems, and solutions exposed in the article; raise criticisms and concerns about current ideas and approaches; suggest current (open) problems with potential solutions; and provide some ideas for future research directions., 36 pages; 6 Figures

Switch-like motifs are among the basic building blocks of biochemical networks. A common motif that can serve as an ultrasensitive switch consists of two enzymes acting antagonistically on a substrate, one making and the other removing a covalent modification. To work as a switch, such covalent modification cycles must be held out of thermodynamic equilibrium by continuous expenditure of energy. Here, we exploit the linear framework for timescale separation to establish tight bounds on the performance of any covalent-modification switch, in terms of the chemical potential difference driving the cycle. The bounds apply to arbitrary enzyme mechanisms, not just Michaelis-Menten, with arbitrary rate constants, and thereby reflect fundamental physical constraints on covalent switching., 29 pages, 6 figures

Ordinary differential equation (ODE) models are widely used to describe chemical or biological processes. This article considers the estimation and assessment of such models on the basis of time-course data. Due to experimental limitations, time-course data are often noisy and some components of the system may not be observed. Furthermore, the computational demands of numerical integration have hindered the widespread adoption of time-course analysis using ODEs. To address these challenges, we explore the efficacy of the recently developed MAGI (MAnifold-constrained Gaussian process Inference) method for ODE inference. First, via a range of examples we show that MAGI is capable of inferring the parameters and system trajectories, including unobserved components, with appropriate uncertainty quantification. Second, we illustrate how MAGI can be used to assess and select different ODE models with time-course data based on MAGI's efficient computation of model predictions. Overall, we believe MAGI is a useful method for the analysis of time-course data in the context of ODE models, which bypasses the need for any numerical integration., Comment: 26 pages, 8 figures, with code supplement

Rapid discovery of new reactions and molecules in recent years has been facilitated by the advancements in high throughput screening, accessibility to a much more complex chemical design space, and the development of accurate molecular modeling frameworks. A holistic study of the growing chemistry literature is, therefore, required that focuses on understanding the recent trends and extrapolating them into possible future trajectories. To this end, several network theory-based studies have been reported that use a directed graph representation of chemical reactions. Here, we perform a study based on representing chemical reactions as hypergraphs where the hyperedges represent chemical reactions and nodes represent the participating molecules. We use a standard reactions dataset to construct a hypernetwork and report its statistics such as degree distributions, average path length, assortativity or degree correlations, PageRank centrality, and graph-based clusters (or communities). We also compute each statistic for an equivalent directed graph representation of reactions to draw parallels and highlight differences between the two. To demonstrate the AI applicability of hypergraph reaction representation, we generate dense hypergraph embeddings and use them in the reaction classification problem. We conclude that the hypernetwork representation is flexible, preserves reaction context, and uncovers hidden insights that are otherwise not apparent in a traditional directed graph representation of chemical reactions.

Control of multistable dynamical system has important applications, from physics to biology. Here, we attack this problem from the perspective of local sensitivity analysis. We develop sensitivity rules to control properties of continuous-time multistable dynamics with simple attractors, namely, the relative size and depth of their basins of attraction. Our parameter control signal is computationally cheap and provides counter-intuitive information about the sensitive parameters to be manipulated in an experimental setting., Comment: Gave an in-depth discussion of and added references to the theory in Section III. Slight change of structure in Control Strategies. Slight changes to figures and discussion of example systems

Perfect adaptation is a phenomenon whereby the output variables of a system can maintain certain values despite external disturbances. Robust perfect adaptation (RPA) refers to an adaptation property that does not require fine-tuning of system parameters. RPA plays a vital role for the survival of living systems in unpredictable environments. However, complex interaction patterns in biochemical systems pose a significant challenge in identifying RPA and associated regulatory mechanisms. The goal of this paper is to present a novel approach for finding all RPA properties that are realized for a generic choice of kinetics for general deterministic chemical reaction systems. This is accomplished by proving that an RPA property is represented by a subnetwork with certain topological features. This connection is exploited to show that these structures generate all kinetics-independent RPA properties, allowing us to systematically identify all RPA properties by enumerating these subnetworks. An efficient method is developed for this enumeration, and we provide a computational package for this purpose. We pinpoint the integral feedback controllers that work in concert to realize each RPA property, casting our results into the familiar control-theoretic paradigm of the Internal Model Principle. We further generalize the regulation problem to the multi-output scenario where the target values belong to a manifold of nonzero dimension, and provide a sufficient condition for this. The present work significantly advances our understanding of regulatory mechanisms that lead to RPA in endogenous biochemical systems, and it also provides rational design principles for synthetic controllers. The present results indicate that an RPA property is essentially equivalent to the existence of a "topological invariant", which is an instance of what we call the "Robust Adaptation is Topological"(RAT) principle., 65 pages, 16 figures

Many reaction networks arising in applications are multistationary, that is, they have the capacity for more than one steady state; while some networks exhibit absolute concentration robustness (ACR), which means that some species concentration is the same at all steady states. Both multistationarity and ACR are significant in biological settings, but only recently has attention focused on the possibility for these properties to coexist. Our main result states that such coexistence in at-most-bimolecular networks (which encompass most networks arising in biology) requires at least $3$ species, $5$ complexes, and $3$ reactions. We prove additional bounds on the number of reactions for general networks based on the number of linear conservation laws. Finally, we prove that, outside of a few exceptional cases, ACR is equivalent to non-multistationarity for bimolecular networks that are small (more precisely, one-dimensional or up to two species). Our proofs involve analyses of systems of sparse polynomials, and we also use classical results from chemical reaction network theory.

There is an unmet need to develop medications or drug combinations which can stop advancement of NASH to liver cirrhosis and HCC. Therefore, identifying key biomarkers based on overall survival and the exploring biological processes involved in the pathogenesis and progression of NASH toward cirrhosis and HCC to improve therapeutic interventions is necessary. The microarray dataset from the GPL13667 platform was downloaded from the Gene Expression Omnibus (GEO). The inclusion criteria for the DEGs included an adjusted p-value 1. In the first step, protein-protein interaction (PPI) network of differentially expressed genes (DEGs) in every three groups (NASH, Cirrhosis and HCC) was constructed using STRING online database. In the second step, the DEGs of each group were imported to Cytoscape software separately, and the genes with Degree >3 were selected. In the third step, the genes with Degree >3 were imported to Gephi software (0.9.2 version) and the genes with betweenness centrality >0 were selected. Venn diagram of these genes was depicted for the NASH, cirrhotic and HCC groups. According to venn diagram, 96 (NASH), 30 (cirrhotic) and 213 (HCC) genes were specifically upregulated (based on inclusion criteria). Among specific genes, 22 (NASH), 5 (cirrhotic) and 82 (HCC) genes were with poor overall survival. The overlap among the 3 groups (NASH, Cirrhosis and HCC) contained 4 upregulated genes HLA-F, HLA-DPA1, TPM1 and YWHAZ. From 4 upregulated genes, only YWHAZ gene was with poor overall survival. The present study detected new candidate genes and key biological processes in NASH, cirrhosis and HCC based on overall survival and using in silico analysis. Therefore, performing in vitro and in vivo researches to verify the results is necessary., 20 pages, 5 figures, 1 tables

The BReast CAncer type 2 susceptibility protein (BRCA2) responds to DNA damage by participating in homology-directed repair. BRCA2 deficiency culminates in defective DNA damage repair (DDR) that when prolonged leads to the accumulation of DNA damage causing cancer or apoptosis. Oxidative stress promotes DNA damage and apoptosis and is a common mechanism through which cardiovascular risk factors lead to endothelial dysfunction and atherosclerosis. Herein, we show that endothelial BRCA2 plays a protective role against atherosclerosis under hypercholesterolemic stress. We successfully generated and characterized endothelial cell (EC)-specific BRCA2 knockout (BRCA2endo) mice. To study the effect of EC-specific BRCA2-loss in atherosclerosis, we generated and characterized BRCA2endo mice on apolipoprotein E null background (ApoE-/-), fed them with high-fat diet (HFD) and evaluated atherosclerosis. Baseline phenotyping of BRCA2endo mice did not show any adverse effects in terms of DNA damage and apoptosis as well as cardiac and metabolic function. However, using HFD-fed apolipoprotein E knockout (ApoE-/-) background, we demonstrated that EC-specific loss of BRCA2 resulted in aortic plaque deposition and splenomegaly. Comparison of RNA sequencing data from aortas of EC-specific BRCA2-deficient ApoE-/- and BRCA2-intact ApoE-/- mice revealed a total of 530 significantly differentially expressed genes with Protein Folding Response and Lipid Metabolism as the most affected pathways. This study provides foundational knowledge regarding BRCA2 status and function in the cardiovascular system, and highlights the potential of BRCA2 as a novel therapeutic target in prevention and treatment of atherosclerosis. Our data indicate that BRCA2 mutation carriers may be at a previously unrecognized risk of atherosclerosis in addition to breast and ovarian cancer.

External flows of energy, entropy, and matter can cause sudden transitions in the stability of biological and industrial systems, fundamentally altering their dynamical function. How might we control and design these transitions in chemical reaction networks? Here, we analyze transitions giving rise to complex behavior in random reaction networks subject to external driving forces. In the absence of driving, we characterize the uniqueness of the steady state and identify the percolation of a giant connected component in these networks as the number of reactions increases. When subject to chemical driving (influx and outflux of chemical species), the steady state can undergo bifurcations, leading to multistability or oscillatory dynamics. By quantifying the prevalence of these bifurcations, we show how chemical driving and network sparsity tend to promote the emergence of these complex dynamics and increased rates of entropy production. We show that catalysis also plays an important role in the emergence of complexity, strongly correlating with the prevalence of bifurcations. Our results suggest that coupling a minimal number of chemical signatures with external driving can lead to features present in biochemical processes and abiogenesis., 12 pages, 4 figures

Logical modeling is a powerful tool in biology, offering a system-level understanding of the complex interactions that govern biological processes. A gap that hinders the scalability of logical models is the need to specify the update function of every vertex in the network depending on the status of its predecessors. To address this, we introduce in this paper the concept of strong regulation, where a vertex is only updated to active/inactive if all its predecessors agree in their influences; otherwise, it is set to ambiguous. We explore the interplay between active, inactive, and ambiguous influences in a network. We discuss the existence of phenotype attractors in such networks, where the status of some of the variables is fixed to active/inactive, while the others can have an arbitrary status, including ambiguous.

We propose a novel alternative approach to our previous work (Ben Hammouda et al., 2023) to improve the efficiency of Monte Carlo (MC) estimators for rare event probabilities for stochastic reaction networks (SRNs). In the same spirit of (Ben Hammouda et al., 2023), an efficient path-dependent measure change is derived based on a connection between determining optimal importance sampling (IS) parameters within a class of probability measures and a stochastic optimal control formulation, corresponding to solving a variance minimization problem. In this work, we propose a novel approach to address the encountered curse of dimensionality by mapping the problem to a significantly lower-dimensional space via a Markovian projection (MP) idea. The output of this model reduction technique is a low-dimensional SRN (potentially even one dimensional) that preserves the marginal distribution of the original high-dimensional SRN system. The dynamics of the projected process are obtained by solving a related optimization problem via a discrete $L^2$ regression. By solving the resulting projected Hamilton-Jacobi-Bellman (HJB) equations for the reduced-dimensional SRN, we obtain projected IS parameters, which are then mapped back to the original full-dimensional SRN system, resulting in an efficient IS-MC estimator for rare events probabilities of the full-dimensional SRN. Our analysis and numerical experiments reveal that the proposed MP-HJB-IS approach substantially reduces the MC estimator variance, resulting in a lower computational complexity in the rare event regime than standard MC estimators.

Gene expression patterns are established by cross-regulating target genes that interpret morphogen gradients. However, as development progresses, morphogen activity is reduced, leaving the emergent pattern without stabilizing positional cues. The pattern then can be deteriorated by the intrinsically noisy biochemical processes acting at the cellular level. But remarkably, the established gene expression patterns remain spatially and temporally stable in many biological systems. Here we combine spatial-stochastic simulations with an enhanced sampling method and a recently developed stability theory to address how spatiotemporal integrity of a gene expression pattern is maintained in developing tissue lacking morphogen gradients. Using a minimal embryo model consisting of spatially coupled biochemical reactor volumes, we study a stripe pattern in which weak cross-repression between nearest neighbor domians alternates with strong repression between next-nearest neighbor domains, inspired by the gap gene system in the Drosophila embryo. We find that fine-tuning of the weak repressive interactions to an optimal level can significantly increase temporal stability of the expression patterns, allowing for stable patterns over developmentally relevant times in the absence of morphogen gradients. The numerically determined optimal parameters closely agree with the predictions of the stability theory. By analizing the dynamics of factors characterizing pattern integrity, we trace back the stability enhancement to the emergence of restoring forces, maintaining the pattern in a meta-stable basin. Altogether our results demonstrate that metastable attractors can emerge as a property of stochastic gene expression patterns even without system-wide positional cues, provided that the gene regulatory interactions shaping the pattern are optimally tuned.

Alterations in the human gut bacteriome can be associated with human health issues, such as type-2 diabetes and cardiovascular disease. Both external and internal factors can drive changes in the composition and in the interactions of the human gut bacteriome, impacting negatively on the host cells. In this paper, we focus on the human gut bacteriome metabolism and we propose a two-layer network system to investigate its dynamics. Furthermore, we develop an in-silico simulation model (virtual GB), allowing us to study the impact of the metabolite exchange through molecular communications in the human gut bacteriome network system. Our results show that the regulation of molecular inputs can strongly affect bacterial population growth and create an unbalanced network, as shown by the shift in the node weights based on the molecular signals that are produced. Additionally, we show that the metabolite molecular communication production is greatly affected when directly manipulating the composition of the human gut bacteriome network in the virtual GB. These results indicate that our human GB interaction model can help to identify hidden behaviors of the human gut bacteriome depending on the molecular signal interactions. Moreover, the virtual GB can support the research and development of novel medical treatments based on the accurate control of bacterial growth and exchange of metabolites.

COVID-19 has created a global pandemic with high morbidity and mortality in 2020. Novel coronavirus (nCoV), also known as Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV2), is responsible for this deadly disease. International Committee on Taxonomy of Viruses (ICTV) has declared that nCoV is highly genetically similar to SARS-CoV epidemic in 2003 (89% similarity). Limited number of clinically validated Human-nCoV protein interaction data is available in the literature. With this hypothesis, the present work focuses on developing a computational model for nCoV-Human protein interaction network, using the experimentally validated SARS-CoV-Human protein interactions. Initially, level-1 and level-2 human spreader proteins are identified in SARS-CoV-Human interaction network, using Susceptible-Infected-Susceptible (SIS) model. These proteins are considered as potential human targets for nCoV bait proteins. A gene-ontology based fuzzy affinity function has been used to construct the nCoV-Human protein interaction network at 99.98% specificity threshold. This also identifies the level-1 human spreaders for COVID-19 in human protein-interaction network. Level-2 human spreaders are subsequently identified using the SIS model. The derived host-pathogen interaction network is finally validated using 7 potential FDA listed drugs for COVID-19 with significant overlap between the known drug target proteins and the identified spreader proteins., Comment: Total number of pages is 12. Total number of figures is 7. Total number of tables is 1

We study the multistationarity for the reaction networks with one-dimensional stoichiometric subspaces, and we focus on the networks admitting finitely many positive steady states. We prove that if a network admits multistationarity, then network has an embedded one-species network with arrow diagram (->,). The inverse is also true if there exist two reactions in the network such that the subnetwork consisting of the two reactions admits at least one and finitely many positive steady states. We also prove that if a network admits at least three positive steady states, then it contains at least three bi-arrow diagrams. More than that, we completely characterize the bi-reaction networks that admit at least three positive steady states., Comment: 31 pages

Recently, molecular relational learning, whose goal is to predict the interaction behavior between molecular pairs, got a surge of interest in molecular sciences due to its wide range of applications. In this work, we propose CMRL that is robust to the distributional shift in molecular relational learning by detecting the core substructure that is causally related to chemical reactions. To do so, we first assume a causal relationship based on the domain knowledge of molecular sciences and construct a structural causal model (SCM) that reveals the relationship between variables. Based on the SCM, we introduce a novel conditional intervention framework whose intervention is conditioned on the paired molecule. With the conditional intervention framework, our model successfully learns from the causal substructure and alleviates the confounding effect of shortcut substructures that are spuriously correlated to chemical reactions. Extensive experiments on various tasks with real-world and synthetic datasets demonstrate the superiority of CMRL over state-of-the-art baseline models. Our code is available at https://github.com/Namkyeong/CMRL., KDD 2023

Given a set of elements, the reactions between them (chemical or otherwise), and certain elements catalysing certain reactions, a Reflexively Autocatalytic F-generated (RAF) set is a subset $R'$ of reactions that is self-generating from a given food set, and with each reaction in $R'$ being catalysed from within $R'$. RAF theory has been applied to various phenomena in theoretical biology, and a key feature of the approach is that it is possible to efficiently identify and classify RAFs within large systems. This is possible because RAFs can be described as the (nonempty) subsets of the reactions that are the fixed points of an (efficiently computable) interior map that operates on subsets of reactions. Although the main generic results concerning RAFs can be derived using just this property, we show that for systems with at least 12 reactions there are generic results concerning RAFs that cannot be proven using the interior operator property alone., 13 pages, 1 figure

Eigenvalue disparity, also known as timescale separation, permits the systematic reduction of deterministic models of enzyme kinetics. Geometric singular perturbation theory, of which eigenvalue disparity is central, provides a coordinate-free framework for deriving reduced mass action models in the deterministic realm. Moreover, homologous deterministic reductions are often employed in stochastic models to reduce the computational complexity required to simulate reactions with the Gillespie algorithm. Interestingly, several detailed studies indicate that timescale separation does not always guarantee the accuracy of reduced stochastic models. In this work, we examine the roles of timescale separation and coordinate transformations in the reduction of the Linear Noise Approximation (LNA) and, unlike previous studies, we do not require the system to be comprised of distinct fast and slow variables. Instead, we adopt a coordinate-free approach. We demonstrate that eigenvalue disparity does not guarantee the accuracy of the reduced LNA, known as the slow scale LNA (ssLNA). However, the inaccuracy of the ssLNA can often be eliminated with a proper coordinate transformation. For planar systems in separated (standard) form, we prove that the error between the variances of the slow variable generated by the LNA and the ssLNA is $\mathcal{O}(\varepsilon)$. We also address a nilpotent Jacobian scenario and use the blow-up method to construct a reduced equation that is accurate near the singular limit in the deterministic regime. However, this reduction in the stochastic regime is far less accurate, which illustrates that eigenvalue disparity plays a central role in stochastic model reduction., 9 Figures, 38 pages

This paper presents a real-time simulation involving ''protozoan-like'' cells that evolve by natural selection in a physical 2D ecosystem. Selection pressure is exerted via the requirements to collect mass and energy from the surroundings in order to reproduce by cell-division. Cells do not have fixed morphologies from birth; they can use their resources in construction projects that produce functional nodes on their surfaces such as photoreceptors for light sensitivity or flagella for motility. Importantly, these nodes act as modular components that connect to the cell's control system via IO channels, meaning that the evolutionary process can replace one function with another while utilising pre-developed control pathways on the other side of the channel. A notable type of node function is the adhesion receptors that allow cells to bind together into multicellular structures in which individuals can share resource and signal to one another. The control system itself is modelled as an artificial neural network that doubles as a gene regulatory network, thereby permitting the co-evolution of form and function in a single data structure and allowing cell specialisation within multicellular groups., Accepted for publication in the proceedings of the 2023 Conference on Artificial Life

Finely-tuned enzymatic pathways control cellular processes, and their dysregulation can lead to disease. Creating predictive and interpretable models for these pathways is challenging because of the complexity of the pathways and of the cellular and genomic contexts. Here we introduce Elektrum, a deep learning framework which addresses these challenges with data-driven and biophysically interpretable models for determining the kinetics of biochemical systems. First, it uses in vitro kinetic assays to rapidly hypothesize an ensemble of high-quality Kinetically Interpretable Neural Networks (KINNs) that predict reaction rates. It then employs a novel transfer learning step, where the KINNs are inserted as intermediary layers into deeper convolutional neural networks, fine-tuning the predictions for reaction-dependent in vivo outcomes. Elektrum makes effective use of the limited, but clean in vitro data and the noisy, yet plentiful in vivo data that captures cellular context. We apply Elektrum to predict CRISPR-Cas9 off-target editing probabilities and demonstrate that Elektrum achieves state-of-the-art performance, regularizes neural network architectures, and maintains physical interpretability., 23 pages, 4 figures

Protein-protein interactions (PPIs) are crucial in various biological processes and their study has significant implications for drug development and disease diagnosis. Existing deep learning methods suffer from significant performance degradation under complex real-world scenarios due to various factors, e.g., label scarcity and domain shift. In this paper, we propose a self-ensembling multigraph neural network (SemiGNN-PPI) that can effectively predict PPIs while being both efficient and generalizable. In SemiGNN-PPI, we not only model the protein correlations but explore the label dependencies by constructing and processing multiple graphs from the perspectives of both features and labels in the graph learning process. We further marry GNN with Mean Teacher to effectively leverage unlabeled graph-structured PPI data for self-ensemble graph learning. We also design multiple graph consistency constraints to align the student and teacher graphs in the feature embedding space, enabling the student model to better learn from the teacher model by incorporating more relationships. Extensive experiments on PPI datasets of different scales with different evaluation settings demonstrate that SemiGNN-PPI outperforms state-of-the-art PPI prediction methods, particularly in challenging scenarios such as training with limited annotations and testing on unseen data., Accepted by IJCAI 2023

RNA structure and functional dynamics play fundamental roles in controlling biological systems. Molecular dynamics simulation, which can characterize interactions at an atomistic level, can advance the understanding on new drug discovery, manufacturing, and delivery mechanisms. However, it is computationally unattainable to support the development of a digital twin for enzymatic reaction network mechanism learning, and end-to-end bioprocess design and control. Thus, we create a hybrid ("mechanistic + machine learning") model characterizing the interdependence of RNA structure and functional dynamics from atomistic to macroscopic levels. To assess the proposed modeling strategy, in this paper, we consider RNA degradation which is a critical process in cellular biology that affects gene expression. The empirical study on RNA lifetime prediction demonstrates the promising performance of the proposed multi-scale bioprocess hybrid modeling strategy., 12 pages, 5 figures

An automata network with $n$ components over a finite alphabet $Q$ of size $q$ is a discrete dynamical system described by the successive iterations of a function $f:Q^n\to Q^n$. In most applications, the main parameter is the interaction graph of $f$: the digraph with vertex set $[n]$ that contains an arc from $j$ to $i$ if $f_i$ depends on input $j$. What can be said on the set $\mathbb{G}(f)$ of the interaction graphs of the automata networks isomorphic to $f$? It seems that this simple question has never been studied. Here, we report some basic facts. First, we prove that if $n\geq 5$ or $q\geq 3$ and $f$ is neither the identity nor constant, then $\mathbb{G}(f)$ always contains the complete digraph $K_n$, with $n^2$ arcs. Then, we prove that $\mathbb{G}(f)$ always contains a digraph whose minimum in-degree is bounded as a function of $q$. Hence, if $n$ is large with respect to $q$, then $\mathbb{G}(f)$ cannot only contain $K_n$. However, we prove that $\mathbb{G}(f)$ can contain only dense digraphs, with at least $\lfloor n^2/4 \rfloor$ arcs., Comment: 20 pages

The rapidly expanding market for regenerative medicines and cell therapies highlights the need to advance the understanding of cellular metabolisms and improve the prediction of cultivation production process for human induced pluripotent stem cells (iPSCs). In this paper, a metabolic kinetic model was developed to characterize underlying mechanisms of iPSC culture process, which can predict cell response to environmental perturbation and support process control. This model focuses on the central carbon metabolic network, including glycolysis, pentose phosphate pathway (PPP), tricarboxylic acid (TCA) cycle, and amino acid metabolism, which plays a crucial role to support iPSC proliferation. Heterogeneous measures of extracellular metabolites and multiple isotopic tracers collected under multiple conditions were used to learn metabolic regulatory mechanisms. Systematic cross-validation confirmed the model's performance in terms of providing reliable predictions on cellular metabolism and culture process dynamics under various culture conditions. Thus, the developed mechanistic kinetic model can support process control strategies to strategically select optimal cell culture conditions at different times, ensure cell product functionality, and facilitate large-scale manufacturing of regenerative medicines and cell therapies., 22 pages, 16 figures

Molecular relational learning, whose goal is to learn the interaction behavior between molecular pairs, got a surge of interest in molecular sciences due to its wide range of applications. Recently, graph neural networks have recently shown great success in molecular relational learning by modeling a molecule as a graph structure, and considering atom-level interactions between two molecules. Despite their success, existing molecular relational learning methods tend to overlook the nature of chemistry, i.e., a chemical compound is composed of multiple substructures such as functional groups that cause distinctive chemical reactions. In this work, we propose a novel relational learning framework, called CGIB, that predicts the interaction behavior between a pair of graphs by detecting core subgraphs therein. The main idea is, given a pair of graphs, to find a subgraph from a graph that contains the minimal sufficient information regarding the task at hand conditioned on the paired graph based on the principle of conditional graph information bottleneck. We argue that our proposed method mimics the nature of chemical reactions, i.e., the core substructure of a molecule varies depending on which other molecule it interacts with. Extensive experiments on various tasks with real-world datasets demonstrate the superiority of CGIB over state-of-the-art baselines. Our code is available at https://github.com/Namkyeong/CGIB., ICML 2023

We consider a generic class of gene circuits affected by nonlinear extrinsic noise. To address this nonlinearity we introduce a general perturbative methodology based on assuming timescale separation between noise and genes dynamics, with fluctuations exhibiting a large but finite correlation time. We apply this methodology to the case of the toggle switch, and by considering biologically relevant log-normal fluctuations, we find that the system exhibits noise-induced transitions. The system becomes bimodal in regions of the parameter space where it would be deterministically monostable. We show that by including higher order corrections our methodology allows one to obtain correct predictions for the occurrence of transitions even for not so large correlation time of the fluctuations, overcoming thereby limitations of previous theoretical approaches. Interestingly we find that at intermediate noise intensities the noise-induced transition in the toggle switch affects one of the genes involved, but not the other one., Comment: 11 pages, 5 figures

The catalytic reaction system (CRS) formalism by Hordijk and Steel is a versatile method to model autocatalytic biochemical reaction networks. It is particularly suited, and has been widely used, to study self-sustainment and self-generation properties. Its distinguishing feature is the explicit assignment of a catalytic function to chemicals that are part of the system. In this work, it is shown that the subsequent and simultaneous catalytic functions gives rise to an algebraic structure of a semigroup with the additional compatible operation of idempotent addition and a partial order. The aim of this article is to demonstrate that such semigroup models are a natural setup to describe and analyze self-sustaining CRS. The basic algebraic properties of the models are established and the notion of the function of any set of chemicals on the whole CRS is made precise. This leads to a natural discrete dynamics on the network, which results from iteratively considering the self-action on a set of chemicals by its own function. The fixed points of this dynamics are proven to correspond to self-sustaining sets of chemicals, which are functionally closed. Finally, as the main application, a theorem on the maximal self-sustaining set and a structure theorem on the lattice of functionally closed self-sustaining sets of chemicals are proven., Comment: 28 pages, 5 figures

When a biological system robustly corrects component-level errors, the direct pressure on component performance declines. Components may become less reliable, maintain more genetic variability, or drift neutrally in design, creating the basis for new forms of organismal complexity. This article links the protection-decay dynamic to other aspects of robust and complex systems. Examples include the hourglass pattern of biological development and Doyle's hourglass architecture for robustly complex systems in engineering. The deeply and densely connected wiring architecture in biology's cellular controls and in machine learning's computational neural networks provide another link. By unifying these seemingly different aspects into a unified framework, we gain a new perspective on robust and complex systems.

The circadian clock can adapt itself to external cues, but the molecular mechanisms and regulatory networks governing circadian oscillations' transient adjustments are still largely unknown. Here we consider the specific case of circadian oscillations transiently responding to a temperature change. Using a framework motivated by Floquet theory, we model the mRNA expression level of the fat body from Drosophila melanogaster following a step change from 25C to 18C. Using the method we infer the adaptation rates of individual genes as they adapt to the new temperature. To deal with heteroskedastic noise and outliers present in the expression data we employ quantile regression and wild bootstrap for significance testing. Model selection with finite-sample corrected Akaike Information Criterion (AICc) is performed additionally for robust inference. We identify several genes with fast transition rates as potential sources of temperature-mediated responses in the circadian system of fruit flies, and the constructed network suggests that the proteasome may play important roles in governing these responses.

Accurately inferring Gene Regulatory Networks (GRNs) is a critical and challenging task in biology. GRNs model the activatory and inhibitory interactions between genes and are inherently causal in nature. To accurately identify GRNs, perturbational data is required. However, most GRN discovery methods only operate on observational data. Recent advances in neural network-based causal discovery methods have significantly improved causal discovery, including handling interventional data, improvements in performance and scalability. However, applying state-of-the-art (SOTA) causal discovery methods in biology poses challenges, such as noisy data and a large number of samples. Thus, adapting the causal discovery methods is necessary to handle these challenges. In this paper, we introduce DiscoGen, a neural network-based GRN discovery method that can denoise gene expression measurements and handle interventional data. We demonstrate that our model outperforms SOTA neural network-based causal discovery methods.

Prokaryotic gene expression is dynamic and noisy, and can lead to phenotypic variations. Previous work has shown that the selection of such variation can be modeled through basic reaction networks described by O.D.E. displaying a bistable behavior. While previous mathematical studies have shown that mono- or bistable behavior depends on the rate of the reactions in the system, no analytical solution of the curves delineating the actual parameter conditions that result in mono or bistability has so far been provided. In this work we provide the first explicit analytical solution for the boundary curve that separates the parameter space defining domains where double positive and double negative feedback loops become bistable.

Causal discovery of genome-scale networks is important for identifying pathways from genes to observable traits - e.g. differences in cell function, disease, drug resistance and others. Causal learners based on graphical models rely on interventional samples to orient edges in the network. However, these models have not been shown to scale up the size of the genome, which are on the order of 1e3-1e4 genes. We introduce a new learner, SP-GIES, that jointly learns from interventional and observational datasets and achieves almost 4x speedup against an existing learner for 1,000 node networks. SP-GIES achieves an AUC-PR score of 0.91 on 1,000 node networks, and scales up to 2,000 node networks - this is 4x larger than existing works. We also show how SP-GIES improves downstream optimal experimental design strategies for selecting interventional experiments to perform on the system. This is an important step forward in realizing causal discovery at scale via autonomous experimental design., To be published in Platform for Advanced Scientific Computing 2023 (PASC23) conference proceedings

We study the generic dimension of the solution set over C^*, R^* and R_{>0} of parametric polynomial systems that consist of linear combinations of monomials scaled by free parameters. We establish a relation between this dimension, Zariski denseness of the set of parameters for which the system has solutions, and the existence of nondegenerate solutions, which enables fast dimension computations. Systems of this form are used to describe the steady states of reaction networks modeled with mass-action kinetics, and as a corollary of our results, we prove that weakly reversible networks have finitely many steady states for generic reaction rate constants and total concentrations.

Crocus sativus is a triploid sterile plant with red stigmas belonging to family of Iridaceae, and sub family Crocoideae. Crocin, picrocrocin, and safranal are three major carotenoid derivatives that are responsible for color, taste and specific aroma of Crocus. Saffron flowers are harvested manually and used as spice, dye or medicinal applications. The natural propagation rate of most geophytes including saffron is relatively low. An in vitro multiplication technique like micropropagation has been used for the propagation of saffron. To understand the efficiency of this alternative and study the molecular basis of apocarotenoid biosynthesis/accumulation, the RT-PCR method was performed on perianth explants that were cultured on MS medium to observe the level of expression of zeaxanthin cleavage dioxygenase (CsZCD) gene during stigma development, and also the impact of soil EC on its expression. The present study was conducted at Plant molecular and physiology Lab, Alzahra University, Tehran, Iran during 2011-2013. Stigma-like structures (SLSs) on calli were collected from immature perianth explants from floral buds of corms that were collected from Ghaen city, and compared to (Torbat-e Haidariye, Mardabad and Shahroud cities) for investigating the impact of different soil EC on CsZCD expression. The results indicated that CsZCD gene was highly expressed in fully developed red SLSs in perianth of cultured samples of Shahroud with the highest salinity. In this research, a close relationship between soil EC and second metabolites regulation is studied. Overall, these results will pave the way for understanding the molecular basis of apocarotenoid biosynthesis and other aspects of stigma development in C. sativus.

Here we introduce simple structures for the analysis of complex hypergraphs, hypergraph animals. These structures are designed to describe the local node neighbourhoods of nodes in hypergraphs. We establish their relationships to lattice animals and network motifs and develop their combinatorial properties, for sparse and uncorrelated hypergraphs. Here we can make use of the tight link of hypergraph animals to partition numbers, which opens up a vast mathematical framework for the analysis of hypergraph animals. We then study their abundances in random hypergraphs. Two transferable insights result from this analysis: (i) it establishes the importance of high-cardinality edges in ensembles of random hypergraphs that are inspired by the classical Erdos-Reny\'i random graphs; and (ii) there is a close connection between degree and hyperedge cardinality in random hypergraphs that shapes the animal abundances and spectra profoundly. Together these two findings imply that we need to spend more effort on investigating and developing suitable conditional ensembles of random hypergraph that can capture real-world structures and their complex dependency structures., Comment: 11 pages, 4 figures

The physical interpretation of the functioning of the adaptive immune system, which has been thoroughly characterized on genetic and molecular levels, provides a unique opportunity to define an adaptive self-organizing biological system in its entirety. This paper describes a configuration space model of immune function, where directed chemical potentials of the system constitute a space of interactions. In the physical sense, the humoral adaptive immune system adjusts the chemical potential of all available antigenic molecules by tuning the chemical potential and organizing the network hierarchy of its sensor-effector molecules, antibodies. The coupling of sensors and effectors allows the system to adjust the thermodynamic activity of antigens and antibodies, while network organization helps minimize chemical potentials and maximize diversity. Mathematically the system couples the variance of Gaussian distributed interaction energies in its interaction space to the exponentially distributed chemical potentials of its effector molecules to maintain its stationary state. This process creates a scale-free network in interaction space, where absolute thermodynamic activity corresponds to node degree. In the thermodynamic interpretation, the system is an ensemble carrying out {mu}N work, adjusting chemical potentials according to the changes in the chemical potentials of the surroundings. The validity of the model is supported by identifying an interaction flexibility index, the corresponding variables in thermodynamics and network science, and by confirming its applicability to the humoral immune system. Overall, this statistical thermodynamics model of adaptive immunity describes how adaptive biological self-organization arises from the maintenance of a scale-free, directed interaction network with fractal topology., 22 pages, 8 figures

Given the large size and complexity of most biochemical regulation and signaling networks, there is a non-trivial relationship between the micro-level logic of component interactions and the observed macro-dynamics. Here we address this issue by formalizing the existing concept of pathway modules, which are sequences of state updates that are guaranteed to occur (barring outside interference) in the dynamics of automata networks after the perturbation of a subset of driver nodes. We present a novel algorithm to automatically extract pathway modules from networks and we characterize the interactions that may take place between modules. This methodology uses only the causal logic of individual node variables (micro-dynamics) without the need to compute the dynamical landscape of the networks (macro-dynamics). Specifically, we identify complex modules, which maximize pathway length and require synergy between their components. This allows us to propose a new take on dynamical modularity that partitions complex networks into causal pathways of variables that are guaranteed to transition to specific states given a perturbation to a set of driver nodes. Thus, the same node variable can take part in distinct modules depending on the state it takes. Our measure of dynamical modularity of a network is then inversely proportional to the overlap among complex modules and maximal when complex modules are completely decouplable from one another in the network dynamics. We estimate dynamical modularity for several genetic regulatory networks, including the Drosophila melanogaster segment-polarity network. We discuss how identifying complex modules and the dynamical modularity portrait of networks explains the macro-dynamics of biological networks, such as uncovering the (more or less) decouplable building blocks of emergent computation (or collective behavior) in biochemical regulation and signaling., 42 pages, 7 figures; updated author information

Developing a mathematical understanding of autocatalysis in chemical reaction networks has both theoretical and practical implications. For a class of autocatalysis, which we term 'stoichiometric autocatalysis', we show that it is possible to classify them in equivalence classes and develop mathematical results about their behavior. We also provide a linear-programming algorithm to exhaustively enumerate them and a scheme to visualize their polyhedral geometry and combinatorics. We then define cluster chemical reaction networks, a framework for coarse-graining realistic chemical reactions using conservation laws. We find that the list of minimal autocatalytic subnetworks in a maximally connected cluster chemical reaction network with one conservation law grows exponentially in the number of species. We end our discussion with open questions concerning autocatalysis and multidisciplinary opportunities for future investigation., 28 pages, 16 figures, 7 tables

Gene regulatory relationships can be abstracted as a gene regulatory network (GRN), which plays a key role in characterizing complex cellular processes and pathways. Recently, graph neural networks (GNNs), as a class of deep learning models, have emerged as a useful tool to infer gene regulatory relationships from gene expression data. However, deep learning models have been found to be vulnerable to noise, which greatly hinders the adoption of deep learning in constructing GRNs, because high noise is often unavoidable in the process of gene expression measurement. Can we preferably prototype a robust GNN for constructing GRNs? In this paper, we give a positive answer by proposing a Quadratic Graph Attention Network (Q-GAT) with a dual attention mechanism. We study the changes in the predictive accuracy of Q-GAT and 9 state-of-the-art baselines by introducing different levels of adversarial perturbations. Experiments in the E. coli and S. cerevisiae datasets suggest that Q-GAT outperforms the state-of-the-art models in robustness. Lastly, we dissect why Q-GAT is robust through the signal-to-noise ratio (SNR) and interpretability analyses. The former informs that nonlinear aggregation of quadratic neurons can amplify useful signals and suppress unwanted noise, thereby facilitating robustness, while the latter reveals that Q-GAT can leverage more features in prediction thanks to the dual attention mechanism, which endows Q-GAT with the ability to confront adversarial perturbation. We have shared our code in https://github.com/Minorway/Q-GAT_for_Robust_Construction_of_GRN for readers' evaluation.

The complexity of mathematical models in biology has rendered model reduction an essential tool in the quantitative biologist’s toolkit. For stochastic reaction networks described using the Chemical Master Equation, commonly used methods include time-scale separation, Linear Mapping Approximation, and state-space lumping. Despite the success of these techniques, they appear to be rather disparate, and at present, no general-purpose approach to model reduction for stochastic reaction networks is known. In this paper, we show that most common model reduction approaches for the Chemical Master Equation can be seen as minimizing a well-known information-theoretic quantity between the full model and its reduction, the Kullback–Leibler divergence defined on the space of trajectories. This allows us to recast the task of model reduction as a variational problem that can be tackled using standard numerical optimization approaches. In addition, we derive general expressions for propensities of a reduced system that generalize those found using classical methods. We show that the Kullback–Leibler divergence is a useful metric to assess model discrepancy and to compare different model reduction techniques using three examples from the literature: an autoregulatory feedback loop, the Michaelis–Menten enzyme system, and a genetic oscillator.

Previous studies have inferred robust stability of reaction networks by utilizing linear programs or iterative algorithms. Such algorithms become tedious or computationally infeasible for large networks. In addition, they operate like black boxes without offering intuition for the structures that are necessary to maintain stability. In this work, we provide several graphical criteria for constructing robust stability certificates, checking robust non-degeneracy, verifying persistence, and establishing global stability. By characterizing a set of stability-preserving graph modifications that includes the enzymatic modification motif, we show that the stability of arbitrarily large nonlinear networks can be examined by simple visual inspection. We show applications of this technique to ubiquitous motifs in systems biology such as post-translational modification (PTM) cycles, the ribosome flow model (RFM), T-cell kinetic proofreading, and others. The results of this paper are dedicated in honor of Eduardo D. Sontag’s seventieth birthday and his pioneering work in nonlinear dynamical systems and mathematical systems biology.

Target Identification by Enzymes (TIE) problem aims to identify the set of enzymes in a given metabolic network, such that their inhibition eliminates a given set of target compounds associated with a disease while incurring minimum damage to the rest of the compounds. This is an NP-complete problem, and thus optimal solutions using classical computers fail to scale to large metabolic networks. In this paper, we consider the TIE problem for identifying drug targets in metabolic networks. We develop the first quantum optimization solution, called QuTIE (Quantum optimization for Target Identification by Enzymes), to this NP-complete problem. We do that by developing an equivalent formulation of the TIE problem in Quadratic Unconstrained Binary Optimization (QUBO) form, then mapping it to a logical graph, which is then embedded on a hardware graph on a quantum computer. Our experimental results on 27 metabolic networks from Escherichia coli, Homo sapiens, and Mus musculus show that QuTIE yields solutions which are optimal or almost optimal. Our experiments also demonstrate that QuTIE can successfully identify enzyme targets already verified in wet-lab experiments for 14 major disease classes.

Many biological processes discriminate between correct and incorrect substrates through the kinetic proofreading mechanism which enables lower error at the cost of higher energy dissipation. Elucidating physicochemical constraints for global minimization of dissipation and error is important for understanding enzyme evolution. Here, we identify theoretically a fundamental error-cost bound which tightly constrains the performance of proofreading networks under any parameter variations preserving the rate discrimination between substrates. The bound is kinetically controlled, i.e. completely determined by the difference between the transition state energies on the underlying free energy landscape. The importance of the bound is analyzed for three biological processes. DNA replication by T7 DNA polymerase is shown to be nearly optimized, i.e. its kinetic parameters place it in the immediate proximity of the error-cost bound. The isoleucyl-tRNA synthetase (IleRS) of E. coli also operates close to the bound, but further optimization is prevented by the need for reaction speed. In contrast, E. coli ribosome operates in a high-dissipation regime, potentially in order to speed up protein production. Together, these findings establish a fundamental error-dissipation relation in biological proofreading networks and provide a theoretical framework for studying error-dissipation trade-off in other systems with biological discrimination., 14 (main text)+ 19 (supporting information) pages, 6+4 figures

Autocatalytic sets are self-catalyzing and self-sustaining chemical reaction networks that are believed to have played an important role in the origin of life. They have been studied extensively both theoretically as well as experimentally. This short note provides (1) a complete and formal definition of autocatalytic sets (or RAF sets), and (2) an efficient algorithm to detect such sets in arbitrary reaction networks. Although both have been presented in various forms in earlier publications, this note serves as a concise and convenient reference., 7 pages, 1 figure

We classify all additive invariants of open Petri nets: these are $\mathbb{N}$-valued invariants which are additive with respect to sequential and parallel composition of open Petri nets. In particular, we prove two classification theorems: one for open Petri nets and one for monically open Petri nets (i.e. open Petri nets whose interfaces are specified by monic maps). Our results can be summarized as follows. The additive invariants of open Petri nets are completely determined by their values on a particular class of single-transition Petri nets. However, for monically open Petri nets, the additive invariants are determined by their values on transitionless Petri nets and all single-transition Petri nets. Our results confirm a conjecture of John Baez (stated during the AMS' 2022 Mathematical Research Communities workshop)., 20 pages

Spatially annotated single-cell datasets provide unprecedented opportunities to dissect cell-cell communication in development and disease. Heterotypic signaling includes interactions between different cell types and is well established in tissue development and spatial organization. Epithelial organization requires several different programs that are tightly regulated. Planar cell polarity is the organization of epithelial cells along the planar axis orthogonal to the apical-basal axis. In this study, we investigate planar cell polarity factors and explore the implications of developmental regulators as malignant drivers. Utilizing cancer systems biology analysis, we derive gene expression network for WNT-ligands (WNT) and their cognate frizzled (FZD) receptors in skin cutaneous melanoma. The profiles supported by unsupervised clustering of multiple-sequence alignments identify ligand-independent signaling and implications for metastatic progression based on the underpinning developmental spatial program. Omics studies and spatial biology connect developmental programs with oncological events and explain key spatial features of metastatic aggressiveness. Dysregulation of prominent planar cell polarity factors such specific representative of the WNT and FZD families in malignant melanoma recapitulates the development program of normal melanocytes but in an uncontrolled and disorganized fashion., Comment: 7 pages, 2 figures

Genetic circuit design is a well-studied problem in synthetic biology. Ever since the first genetic circuits─the repressilator and the toggle switch─were designed and implemented, many advances have been made in this area of research. The current review systematically organizes a number of key works in this domain by employing the versatile framework of generalized morphological analysis. Literature in the area has been mapped on the basis of (a) the design methodologies used, ranging from brute-force searches to control-theoretic approaches, (b) the modeling techniques employed, (c) various circuit functionalities implemented, (d) key design characteristics, and (e) the strategies used for the robust design of genetic circuits. We conclude our review with an outlook on multiple exciting areas for future research, based on the systematic assessment of key research gaps that have been readily unravelled by our analysis framework.

Mass spectrometry (MS) based omics data analysis require significant time and resources. To date, few parallel algorithms have been proposed for deducing peptides from mass spectrometry-based data. However, these parallel algorithms were designed, and developed when the amount of data that needed to be processed was smaller in scale. In this paper, we prove that the communication bound that is reached by the \emph{existing} parallel algorithms is $\Omega(mn+2r\frac{q}{p})$, where $m$ and $n$ are the dimensions of the theoretical database matrix, $q$ and $r$ are dimensions of spectra, and $p$ is the number of processors. We further prove that communication-optimal strategy with fast-memory $\sqrt{M} = mn + \frac{2qr}{p}$ can achieve $\Omega({\frac{2mnq}{p}})$ but is not achieved by any existing parallel proteomics algorithms till date. To validate our claim, we performed a meta-analysis of published parallel algorithms, and their performance results. We show that sub-optimal speedups with increasing number of processors is a direct consequence of not achieving the communication lower-bounds. We further validate our claim by performing experiments which demonstrate the communication bounds that are proved in this paper. Consequently, we assert that next-generation of \emph{provable}, and demonstrated superior parallel algorithms are urgently needed for MS based large systems-biology studies especially for meta-proteomics, proteogenomic, microbiome, and proteomics for non-model organisms. Our hope is that this paper will excite the parallel computing community to further investigate parallel algorithms for highly influential MS based omics problems., Comment: 20 pages, 6 figures, preprint

Genealogical networks (i.e. family trees) are of growing interest, with the largest known data sets now including well over one billion individuals. Interest in family history also supports an 8.5 billion dollar industry whose size is projected to double within 7 years [FutureWise report HC-1137]. Yet little mathematical attention has been paid to the complex network properties of genealogical networks, especially at large scales. The structure of genealogical networks is of particular interest due to the practice of forming unions, e.g. marriages, that are typically well outside one’s immediate family. In most other networks, including other social networks, no equivalent restriction exists on the distance at which relationships form. To study the effect this has on genealogical networks we use persistent homology to identify and compare the structure of 101 genealogical and 31 other social networks. Specifically, we introduce the notion of a network’s persistence curve, which encodes the network’s set of persistence intervals. We find that the persistence curves of genealogical networks have a distinct structure when compared to other social networks. This difference in structure also extends to subnetworks of genealogical and social networks suggesting that, even with incomplete data, persistent homology can be used to meaningfully analyze genealogical networks. Here we also describe how concepts from genealogical networks, such as common ancestor cycles, are represented using persistent homology. We expect that persistent homology tools will become increasingly important in genealogical exploration as popular interest in ancestry research continues to expand.