Search

Fluctuations are inherent to biological systems, arising from the stochastic nature of molecular interactions, and influence various aspects of system behavior, stability, and robustness. These fluctuations can be categorized as intrinsic, stemming from the system's inherent structure and dynamics, and extrinsic, arising from external factors, such as temperature variations. Understanding the interplay between these fluctuations is crucial for obtaining a comprehensive understanding of biological phenomena. However, studying these effects poses significant computational challenges. In this study, we used an underexplored methodology to analyze the effect of extrinsic fluctuations in stochastic systems using ordinary differential equations instead of solving the Master Equation with stochastic parameters. By incorporating temperature fluctuations into reaction rates, we explored the impact of extrinsic factors on system dynamics. We constructed a master equation and calculated the equations for the dynamics of the first two moments, offering computational efficiency compared with directly solving the chemical master equation. We applied this approach to analyze a biological oscillator, focusing on the p53 model and its response to temperature-induced extrinsic fluctuations. Our findings underscore the impact of extrinsic fluctuations on the nature of oscillations in biological systems, with alterations in oscillatory behavior depending on the characteristics of extrinsic fluctuations. We observed an increased oscillation amplitude and frequency of the p53 concentration cycle. This study provides valuable insights into the effects of extrinsic fluctuations on biological oscillations and highlights the importance of considering them in more complex systems to prevent unwanted scenarios related to health issues., Comment: 21 pages, 9 figures

Models of reaction networks within interacting compartments (RNIC) are a generalization of stochastic reaction networks. It is most natural to think of the interacting compartments as ``cells'' that can appear, degrade, split, and even merge, with each cell containing an evolving copy of the underlying stochastic reaction network. Such models have a number of parameters, including those associated with the internal chemical model and those associated with the compartment interactions, and it is natural to want efficient computational methods for the numerical estimation of sensitivities of model statistics with respect to these parameters. Motivated by the extensive work on computational methods for parametric sensitivity analysis in the context of stochastic reaction networks over the past few decades, we provide a number of methods in the RNIC setting. Provided methods include the (unbiased) Girsanov transformation method (also called the Likelihood Ratio method) and a number of coupling methods for the implementation of finite differences. We provide several numerical examples and conclude that the method associated with the ``Split Coupling'' provides the most efficient algorithm. This finding is in line with the conclusions from the work related to sensitivity analysis of standard stochastic reaction networks. We have made all of the Matlab code used to implement the various methods freely available for download., Comment: 31 pages, a number of images

The biological targets of traditional Chinese medicine (TCM) are the core effectors mediating the interaction between TCM and the human body. Identification of TCM targets is essential to elucidate the chemical basis and mechanisms of TCM for treating diseases. Given the chemical complexity of TCM, both in silico high-throughput drug-target interaction predicting models and biological profile-based methods have been commonly applied for identifying TCM targets based on the structural information of TCM chemical components and biological information, respectively. However, the existing methods lack the integration of TCM chemical and biological information, resulting in difficulty in the systematic discovery of TCM action pathways. To solve this problem, we propose a novel target identification model NP-TCMtarget to explore the TCM target path by combining the overall chemical and biological profiles. First, NP-TCMtarget infers TCM effect targets by calculating associations between drug/disease inducible gene expression profiles and specific gene signatures for 8,233 targets. Then, NP-TCMtarget utilizes a constructed binary classification model to predict binding targets of herbal ingredients. Finally, we can distinguish TCM direct and indirect targets by comparing the effect targets and binding targets to establish the action pathways of herbal components-direct targets-indirect targets by mapping TCM targets in the biological molecular network. We apply NP-TCMtarget to the formula XiaoKeAn to demonstrate the power of revealing the action pathways of herbal formula. We expect that this novel model could provide a systematic framework for exploring the molecular mechanisms of TCM at the target level. NP-TCMtarget is available at http://www.bcxnfz.top/NP-TCMtarget., Comment: 29 pages, 4 figures

The description of gene interactions that constantly occur in the cellular environment is an extremely challenging task due to an immense number of degrees of freedom and incomplete knowledge about microscopic details. Hence, a coarse-grained and rather powerful modeling of such dynamics is provided by Boolean Networks (BNs). BNs are dynamical systems composed of Boolean agents and a record of their possible interactions over time. Stable states in these systems are called attractors which are closely related to the cellular expression of biological phenotypes. Identifying the full set of attractors is, therefore, of substantial biological interest. However, for conventional high-performance computing, this problem is plagued by an exponential growth of the dynamic state space. Here, we demonstrate a novel quantum search algorithm inspired by Grover's algorithm to be implemented on quantum computing platforms. The algorithm performs an iterative suppression of states belonging to basins of previously discovered attractors from a uniform superposition, thus increasing the amplitudes of states in basins of yet unknown attractors. This approach guarantees that a new attractor state is measured with each iteration of the algorithm, an optimization not currently achieved by any other algorithm in the literature. Tests of its resistance to noise have also shown promising performance on devices from the current Noise Intermediate Scale Quantum Computing (NISQ) era., Comment: 9 pages, 6 figures

Gene expression prediction plays a vital role in transcriptome-wide association studies (TWAS), which seek to establish associations between tissue gene expression and complex traits. Traditional models rely on genetic variants in close genomic proximity to the gene of interest to predict the genetic component of gene expression. In this study, we propose a novel approach incorporating distal genetic variants acting through gene regulatory networks (GRNs) into gene expression prediction models, in line with the omnigenic model of complex trait inheritance. Using causal and coexpression GRNs reconstructed from genomic and transcriptomic data and modeling the data as a Bayesian network jointly over genetic variants and genes, inference of gene expression from observed genotypic data is achieved through a two-step process. Initially, the expression level of each gene in the network is predicted using its local genetic variants. The residuals, calculated as the differences between the observed and predicted expression levels, are then modeled using the genotype information of parent and/or grandparent nodes in the GRN. The final predicted expression level of the gene is obtained by summing the predictions from the local variants model and the residual model, effectively incorporating both local and distal genetic influences. Using various regularized regression techniques for parameter estimation, we found that GRN-based gene expression prediction outperformed the traditional local-variant approach on simulated data from the DREAM5 Systems Genetics Challenge and real data from the Geuvadis study and an eQTL mapping study in yeast. This study provides important insights into the challenge of gene expression prediction for TWAS. It reaffirms the importance of GRNs for understanding the genetic effects on gene expression and complex traits more generally., Comment: 19 pages, 3 figures, 2 tables + 11 pages, 8 figures, 1 table supplementary information; code is available on GitHub at https://github.com/guutama/GRN-TI.git

Research organisms provide invaluable insights into human biology and diseases, serving as essential tools for functional experiments, disease modeling, and drug testing. However, evolutionary divergence between humans and research organisms hinders effective knowledge transfer across species. Here, we review state-of-the-art methods for computationally transferring knowledge across species, primarily focusing on methods that utilize transcriptome data and/or molecular networks. We introduce the term "agnology" to describe the functional equivalence of molecular components regardless of evolutionary origin, as this concept is becoming pervasive in integrative data-driven models where the role of evolutionary origin can become unclear. Our review addresses four key areas of information and knowledge transfer across species: (1) transferring disease and gene annotation knowledge, (2) identifying agnologous molecular components, (3) inferring equivalent perturbed genes or gene sets, and (4) identifying agnologous cell types. We conclude with an outlook on future directions and several key challenges that remain in cross-species knowledge transfer.

New experimental methods make it possible to measure the expression levels of many genes, simultaneously, in snapshots from thousands or even millions of individual cells. Current approaches to analyze these experiments involve clustering or low-dimensional projections. Here we use the principle of maximum entropy to obtain a probabilistic description that captures the observed presence or absence of mRNAs from hundreds of genes in cells from the mammalian brain. We construct the Ising model compatible with experimental means and pairwise correlations, and validate it by showing that it gives good predictions for higher-order statistics. We notice that the probability distribution of cell states has many local maxima. By labeling cell states according to the associated maximum, we obtain a cell classification that agrees well with previous results that use traditional clustering techniques. Our results provide quantitative descriptions of gene expression statistics and interpretable criteria for defining cell classes, supporting the hypothesis that cell classes emerge from the collective interaction of gene expression levels., Comment: 14 pages, 21 figures

Elucidating gene function is one of the ultimate goals of plant science. Despite this, only ~15% of all genes in the model plant Arabidopsis thaliana have comprehensively experimentally verified functions. While bioinformatical gene function prediction approaches can guide biologists in their experimental efforts, neither the performance of the gene function prediction methods nor the number of experimental characterisation of genes has increased dramatically in recent years. In this review, we will discuss the status quo and the trajectory of gene function elucidation and outline the recent advances in gene function prediction approaches. We will then discuss how recent artificial intelligence advances in large language models and knowledge graphs can be leveraged to accelerate gene function predictions and keep us updated with scientific literature.

We prove that nested canalizing functions are the minimum-sensitivity Boolean functions for any activity ratio and we determine the functional form of this boundary which has a nontrivial fractal structure. We further observe that the majority of the gene regulatory functions found in known biological networks (submitted to the Cell Collective database) lie on the line of minimum sensitivity which paradoxically remains largely in the unstable regime. Our results provide a quantitative basis for the argument that an evolutionary preference for nested canalizing functions in gene regulation (e.g., for higher robustness) and for elasticity of gene activity are sufficient for concentration of such systems near the "edge of chaos." The original structure of gene regulatory networks is unknown due to the undiscovered functions of some genes. Most gene function discovery approaches make use of unsupervised clustering or classification methods that discover and exploit patterns in gene expression profiles. However, existing knowledge in the field derives from multiple and diverse sources. Incorporating this know-how for novel gene function prediction can, therefore, be expected to improve such predictions. We here propose a function-specific novel gene discovery tool that uses a semi-supervised autoencoder. Our method is thus able to address the needs of a modern researcher whose expertise is typically confined to a specific functional domain. Lastly, the dynamics of unorthodox learning approaches like biologically plausible learning algorithms are investigated and found to exhibit a general form of Einstein relation.

Understanding the mechanisms of interactions within cells, tissues, and organisms is crucial to driving developments across biology and medicine. Mathematical modeling is an essential tool for simulating biological systems and revealing biochemical regulatory mechanisms. Building on experiments, mechanistic models are widely used to describe small-scale intracellular networks and uncover biochemical mechanisms in healthy and diseased states. The rapid development of high-throughput sequencing techniques and computational tools has recently enabled models that span multiple scales, often integrating signaling, gene regulatory, and metabolic networks. These multiscale models enable comprehensive investigations of cellular networks and thus reveal previously unknown disease mechanisms and pharmacological interventions. Here, we review systems biology models from classical mechanistic models to larger, multiscale models that integrate multiple layers of cellular networks. We introduce several examples of models of hypertrophic cardiomyopathy, exercise, and cancer cell proliferation. Additionally, we discuss methods that increase the certainty and accuracy of model predictions. Integrating multiscale models has become a powerful tool for understanding disease and inspiring drug discoveries by incorporating omics data within the cell and across tissues and organisms., Comment: 32 pages, 4 figures

Stress generation by the actin cytoskeleton shapes cells and tissues. Despite impressive progress in live imaging and quantitative physical descriptions of cytoskeletal network dynamics, the connection between processes at molecular scales and cell-scale spatio-temporal patterns is still unclear. Here we review studies reporting acto-myosin clusters of micrometer size and with lifetimes of several minutes in a large number of organisms ranging from fission yeast to humans. Such structures have also been found in reconstituted systems in vitro and in theoretical analysis of cytoskeletal dynamics. We propose that tracking these clusters can serve as a simple readout for characterising living matter. Spatio-temporal patterns of clusters could serve as determinants of morphogenetic processes that play similar roles in diverse organisms.

Bacterial second messengers are crucial for transmitting environmental information to cellular responses. However, quantifying their information transmission capacity remains challenging. Here, we engineer an isolated cAMP signaling channel in Pseudomonas aeruginosa using targeted gene knockouts, optogenetics, and a fluorescent cAMP probe. This design allows precise optical control and real-time monitoring of cAMP dynamics. By integrating experimental data with information theory, we reveal an optimal frequency for light-mediated cAMP signaling that maximizes information transmission, reaching about 40 bits/h. This rate correlates strongly with cAMP degradation kinetics and employs a two-state encoding scheme. Our findings suggest a mechanism for fine-tuned regulation of multiple genes through temporal encoding of second messenger signals, providing new insights into bacterial adaptation strategies. This approach offers a framework for quantifying information processing in cellular signaling systems., Comment: 33 pages, 4 figures

Notch signaling is a ubiquitous and versatile intercellular signaling system that drives collective behaviors and pattern formation in biological tissues. During embryonic development, Notch is involved in generation of collective biochemical oscillations that form the vertebrate body segments, and its failure results in embryonic defects. Notch ligands of the Delta family are key components of this collective rhythm, but it is unclear how different Delta ligands with distinct properties contribute to relaying information among cells. Motivated by the zebrafish segmentation clock, in this work we propose a theory describing interactions between biochemical oscillators, where Notch receptor is bound by both oscillatory and nonoscillatory Delta ligands. Based on previous in vitro binding studies, we first consider Notch activation by Delta dimers. This hypothesis is consistent with experimental observations in conditions of perturbed Notch signaling. Then we test an alternative hypothesis where Delta monomers directly bind and activate Notch, and show that this second model can also describe the experimental observations. We show that these two hypotheses assign different roles for a non-oscillatory ligand, as a binding partner or as a baseline signal. Finally, we discuss experiments to distinguish between the two scenarios. Broadly, this work highlights how a multiplicity of ligands may be harnessed by a signaling system to generate versatile responses.

Living systems are usually maintained out of equilibrium and exhibit complex dynamical behaviors. The external energy supply often comes from chemical fluxes that can keep some species concentrations constant. Furthermore, the properties of the underlying chemical reaction networks (CRNs) are also instrumental in establishing robust biological functioning. Hence, capturing the emergent complexity of living systems and the role of their non-equilibrium nature is fundamental to uncover constraints and properties of the CRNs underpinning their functions. In particular, while kinetics plays a key role in shaping detailed dynamical phenomena, the range of operations of any CRN must be fundamentally constrained by thermodynamics, as they necessarily operate with a given energy budget. Here, we derive universal thermodynamic upper and lower bounds for the accessible space of species concentrations in a generic CRN. The resulting region determines the "thermodynamic space" of the CRN, a concept we introduce in this work. Moreover, we obtain similar bounds also for the affinities, shedding light on how global thermodynamic properties can limit local non-equilibrium quantities. We illustrate our results in two paradigmatic examples, the Schl\"ogl model for bistability and a minimal self-assembly process, demonstrating how the onset of complex behaviors is intimately tangled with the presence of non-equilibrium driving. In summary, our work unveils the exact form of the accessible space in which a CRN must work as a function of its energy budget, shedding light on the non-equilibrium origin of a variety of phenomena, from amplification to pattern formation. Ultimately, by providing a general tool for analyzing CRNs, the presented framework constitutes a stepping stone to deepen our ability to predict complex out-of-equilibrium behaviors and design artificial chemical reaction systems.

Markov chains on the non-negative quadrant of dimension $d$ are often used to model the stochastic dynamics of the number of $d$ entities, such as $d$ chemical species in stochastic reaction networks. The infinite state space poses technical challenges, and the boundary of the quadrant can have a dramatic effect on the long term behavior of these Markov chains. For instance, the boundary can slow down the convergence speed of an ergodic Markov chain towards its stationary distribution due to the extinction or the lack of an entity. In this paper, we quantify this slow-down for a class of stochastic reaction networks and for more general Markov chains on the non-negative quadrant. We establish general criteria for such a Markov chain to exhibit a power-law lower bound for its mixing time. The lower bound is of order $|x|^\theta$ for all initial state $x$ on a boundary face of the quadrant, where $\theta$ is characterized by the local behavior of the Markov chain near the boundary of the quadrant. A better understanding of how these lower bounds arise leads to insights into how the structure of chemical reaction networks contributes to slow-mixing., Comment: 39 pages

Homeostasis occurs in a biological system when a chosen output variable remains approximately constant despite changes in an input variable. In this work we specifically focus on biological systems which may be represented as chemical reaction networks and consider their infinitesimal homeostasis, where the derivative of the input-output function is zero. The specific challenge of chemical reaction networks is that they often obey various conservation laws complicating the standard input-output analysis. We derive several results that allow to verify the existence of infinitesimal homeostasis points both in the absence of conservation and under conservation laws where conserved quantities serve as input parameters. In particular, we introduce the notion of infinitesimal concentration robustness, where the output variable remains nearly constant despite fluctuations in the conserved quantities. We provide several examples of chemical networks which illustrate our results both in deterministic and stochastic settings., Comment: 5 figures

Motivation: In systems biology, modelling strategies aim to decode how molecular components interact to generate dynamical behaviour. Boolean modelling is more and more used, but the description of the dynamics from two-levels components may be too limited to capture certain dynamical properties. %However, in Boolean models, the description of the dynamics may be too limited to capture certain dynamical properties. Multivalued logical models can overcome this limitation by allowing more than two levels for each component. However, multivaluing a Boolean model is challenging. Results: We present MRBM, a method for efficiently identifying the components of a Boolean model to be multivalued in order to capture specific fixed-point reachabilities in the asynchronous dynamics. To this goal, we defined a new updating scheme locating reachability properties in the most permissive dynamics. MRBM is supported by mathematical demonstrations and illustrated on a toy model and on two models of stem cell differentiation.

We develop a framework describing the dynamics and thermodynamics of open non-ideal reaction-diffusion systems, which embodies Flory-Huggins theories of mixtures and chemical reaction network theories. Our theory elucidates the mechanisms underpinning the emergence of self-organized dissipative structures in these systems. It evaluates the dissipation needed to sustain and control them, discriminating the contributions from each reaction and diffusion process with spatial resolution. It also reveals the role of the reaction network in powering and shaping these structures. We identify particular classes of networks in which diffusion processes always equilibrate within the structures, while dissipation occurs solely due to chemical reactions. The spatial configurations resulting from these processes can be derived by minimizing a kinetic potential, contrasting with the minimization of the thermodynamic free energy in passive systems. This framework opens the way to investigating the energetic cost of phenomena such as liquid-liquid phase separation, coacervation, and the formation of biomolecular condensates.

Interrogation of biological pathways is an integral part of omics data analysis. Large language models (LLMs) enable the generation of custom pathways and gene sets tailored to specific scientific questions. These targeted sets are significantly smaller than traditional pathway enrichment analysis libraries, reducing multiple hypothesis testing and potentially enhancing statistical power. Lomics (Large Language Models for Omics Studies) v1.0 is a python-based bioinformatics toolkit that streamlines the generation of pathways and gene sets for transcriptomic analysis. It operates in three steps: 1) deriving relevant pathways based on the researcher's scientific question, 2) generating valid gene sets for each pathway, and 3) outputting the results as .GMX files. Lomics also provides explanations for pathway selections. Consistency and accuracy are ensured through iterative processes, JSON format validation, and HUGO Gene Nomenclature Committee (HGNC) gene symbol verification. Lomics serves as a foundation for integrating LLMs into omics research, potentially improving the specificity and efficiency of pathway analysis.

In pharmaceutical research, the strategy of drug repurposing accelerates the development of new therapies while reducing R&D costs. Network pharmacology lays the theoretical groundwork for identifying new drug indications, and deep graph models have become essential for their precision in mapping complex biological networks. Our study introduces an advanced graph model that utilizes graph convolutional networks and tensor decomposition to effectively predict signed chemical-gene interactions. This model demonstrates superior predictive performance, especially in handling the polar relations in biological networks. Our research opens new avenues for drug discovery and repurposing, especially in understanding the mechanism of actions of drugs.

Exchange of molecules allows cells to exchange information. How robust is the information to changes in cell parameters? We use a mapping between the stochastic dynamics of two cells sharing a stimulatory molecule, and parameters akin to an extension of Landau's equilibrium phase transition theory. We show that different single-cell dynamics lead to the same dynamical response -- a flexibility that cells can use. The companion equilibrium Landau model behaves similarly, thereby describing the dynamics of information in a broad class of models with coupled order parameters., Comment: Main text: 6 pages, 3 figures; Supplement: 2 page, 1 figure

How much information does a cell inherit from its ancestors beyond its genetic sequence? What are the epigenetic mechanisms that allow this? Despite the rise in available epigenetic data, how such information is inherited through the cell cycle is still not fully understood. Here, we develop and analyse a simple mathematical model for histone-based epigenetic information that describes how daughter cells can recapitulate the gene expression profiles of their parent. We consider the dynamics of histone modifications during the cell cycle deterministically but also incorporate the largest stochastic element: DNA replication, where histones are randomly distributed between the two daughter DNA strands. This hybrid stochastic-deterministic approach enables an analytic derivation of the switching rate, i.e., the frequency of loss-of-memory events due to replication. While retaining great simplicity, the model can recapitulate experimental switching rate data, establishing its biological importance as a framework to quantitatively study epigenetic inheritance., Comment: 15 pages, 3 figures

Boolean networks constitute relevant mathematical models to study the behaviours of genetic and signalling networks. These networks define regulatory influences between molecular nodes, each being associated to a Boolean variable and a regulatory (local) function specifying its dynamical behaviour depending on its regulators. However, existing data is mostly insufficient to adequately parametrise a model, that is to uniquely define a regulatory function for each node. With the intend to support model parametrisation, this paper presents results on the set of Boolean functions compatible with a given regulatory structure, i.e. the partially ordered set of monotone non-degenerate Boolean functions. More precisely, we present original rules to obtain the direct neighbours of any function of this set. Besides a theoretical interest, presented results will enable the development of more efficient methods for Boolean network synthesis and revision, benefiting from the progressive exploration of the vicinity of regulatory functions., Comment: arXiv admin note: text overlap with arXiv:1901.07623

We present the first use of machine learning tools to predict multistationarity of reaction networks. Chemical Reaction Networks (CRNs) are the mathematical formulation of how the quantities associated to a set of species (molecules, proteins, cells, or animals) vary as time passes with respect to their interactions with each other. Their mathematics does not describe just chemical reactions but many other areas of the life sciences such as ecology, epidemiology, and population dynamics. We say a CRN is at a steady state when the concentration (or number) of species do not vary anymore. Some CRNs do not attain a steady state while some others may have more than one possible steady state. The CRNs in the later group are called multistationary. Multistationarity is an important property, e.g. switch-like behaviour in cells needs multistationarity to occur. Existing algorithms to detect whether a CRN is multistationary or not are either extremely expensive or restricted in the type of CRNs they can be used on, motivating a new machine learning approach. We address the problem of representing variable-length CRN data to machine learning models by developing a new graph representation of CRNs for use with graph learning algorithms. We contribute a large dataset of labelled fully open CRNs whose production necessitated the development of new CRN theory. Then we present experimental results on the training and testing of a graph attention network model on this dataset, showing excellent levels of performance. We finish by testing the model predictions on validation data produced independently, demonstrating generalisability of the model to different types of CRN., Comment: 39 pages, 10 Figures, the dataset and code related to this manuscript is available at the Zenodo link given inside the paper

Proteins can be regarded as thermal nanosensors in an intra-body network. Upon being stimulated by Terahertz (THz) frequencies that match their vibrational modes, protein molecules experience resonant absorption and dissipate their energy as heat, undergoing a thermal process. This paper aims to analyze the effect of THz signaling on the protein heat dissipation mechanism. We therefore deploy a mathematical framework based on the heat diffusion model to characterize how proteins absorb THz-electromagnetic (EM) energy from the stimulating EM fields and subsequently release this energy as heat to their immediate surroundings. We also conduct a parametric study to explain the impact of the signal power, pulse duration, and interparticle distance on the protein thermal analysis. In addition, we demonstrate the relationship between the change in temperature and the opening probability of thermally-gated ion channels. Our results indicate that a controlled temperature change can be achieved in an intra-body environment by exciting protein particles at their resonant frequencies. We further verify our results numerically using COMSOL Multiphysics and introduce an experimental framework that assesses the effects of THz radiation on protein particles. We conclude that under controlled heating, protein molecules can serve as hotspots that impact thermally-gated ion channels. Through the presented work, we infer that the heating process can be engineered on different time and length scales by controlling the THz-EM signal input., Comment: Accepted for publication in the IEEE Transactions on NanoBioscience

Mathematical models of gene regulatory networks are widely used to study cell fate changes and transcriptional regulation. When designing such models, it is important to accurately account for sources of stochasticity. However, doing so can be computationally expensive and analytically untractable, posing limits on the extent of our explorations and on parameter inference. Here, we explore this challenge using the example of a simple auto-negative feedback motif, in which we incorporate stochastic variation due to transcriptional bursting and noise from finite copy numbers. We find that transcriptional bursting may change the qualitative dynamics of the system by inducing oscillations when they would not otherwise be present, or by magnifying existing oscillations. We describe multiple levels of approximation for the model in the form of differential equations, piecewise deterministic processes, and stochastic differential equations. Importantly, we derive how the classical chemical Langevin equation can be extended to include a noise term representing transcriptional bursting. This approximation drastically decreases computation times and allows us to analytically calculate properties of the dynamics, such as their power spectrum. We explore when these approximations break down and provide recommendations for their use. Our analysis illustrates the importance of accounting for transcriptional bursting when simulating gene regulatory network dynamics and provides recommendations to do so with computationally efficient methods.

Complex systems can be advantageously modeled by formal reaction systems (RS), a.k.a. chemical reaction networks in chemistry. Reaction-based models can indeed be interpreted in a hierarchy of semantics, depending on the question at hand, most notably by Ordinary Differential Equations (ODEs), Continuous Time Markov Chains (CTMCs), discrete Petri nets and asynchronous Boolean transition systems. The last three semantics can be easily related in the framework of abstract interpretation. The first two are classically related by Kurtz's limit theorem which states that if reactions are density-dependent families, then, as the volume goes to infinity, the mean reactant concentrations of the CTMC tends towards the solution of the ODE. In the more realistic context of bounded volumes, it is easy to show, by moment closure, that the restriction to reactions with at most one reactant ensures similarly that the mean of the CTMC trajectories is equal to the solution of the ODE at all time points. In this paper, we generalize that result in presence of polyreactant reactions, by introducing the Stoichiometric Influence and Modification Graph (SIMG) of an RS, and by showing that the equality between the two interpretations holds for the variables that belong to distinct SIMG ancestors of polyreactant reactions. We illustrate this approach with several examples. Evaluation on BioModels reveals that the condition for all variables is satisfied on models with no polymolecular reaction only. However, our theorem can be applied selectively to certain variables of the model to provide insights into their behaviour within more complex systems. Interestingly, we also show that the equality holds for a basic oscillatory RS implementing the sine and cosine functions of time.

Inferring parameters of models of biochemical kinetics from single-cell data remains challenging because of the uncertainty arising from the intractability of the likelihood function of stochastic reaction networks. Such uncertainty falls beyond current error quantification measures, which focus on the effects of finite sample size and identifiability but lack theoretical guarantees when likelihood approximations are needed. Here, we propose an inference method for stochastic reaction networks with nonlinear and rational propensities at steady state that provides bounds on the parameters via convex optimisation over sets constrained by moment equations and moment matrices. Our approach takes observations from the stochastic reaction network and forms moment intervals, which are then used to constrain parameters through convex sets. The bounds on the parameters contain the true parameters under the condition that the moment intervals contain the true stationary moments, thus providing uncertainty quantification and error guarantees. Our approach does not need to predict moments and distributions for given parameters (i.e., it avoids solving or simulating the forward problem), and hence circumvents intractable likelihood computations or computationally expensive simulations. We demonstrate its use for uncertainty quantification, data integration and prediction of latent species statistics through synthetic data from common nonlinear biochemical models including the Schl\"ogl model, the toggle switch and post-transcriptional regulation.

Metabolic models condense biochemical knowledge about organisms in a structured and standardised way. As large-scale network reconstructions are readily available for many organisms of interest, genome-scale models are being widely used among modellers and engineers. However, these large models can be difficult to analyse and visualise, and occasionally generate hard-to-interpret or even biologically unrealistic predictions. Out of the thousands of enzymatic reactions in a typical bacterial metabolism, only a few hundred comprise the metabolic pathways essential to produce energy carriers and biosynthetic precursors. These pathways carry relatively high flux, are central to maintaining and reproducing the cell, and provide precursors and energy to engineered metabolic pathways. Here, focusing on these central metabolic subsystems, we present a manually-curated medium-scale model of energy and biosynthesis metabolism for the well-studied prokaryote Escherichia coli K-12 MG1655. The model is a sub-network of the most recent genome-scale reconstruction, iML1515, and comes with an updated layer of database annotations, as well as a range of metabolic maps for visualisation. We enriched the stoichiometric network with extensive biological information and quantitative data, enhancing the scope and applicability of the model. In addition, here we assess the properties of this model in relation to its genome-scale parent and demonstrate the use of the network and supporting data in various scenarios, including enzyme-constrained flux balance analysis, elementary flux mode analysis, and thermodynamic analysis. Overall, we believe this model holds the potential to become a reference medium-scale metabolic model for E. coli.

Bacterial sensor systems can be used for the detection and measurement of molecular signal concentrations. The dynamics of the sensor directly depend on the biological properties of the bacterial sensor cells; manipulation of these features in the wet lab enables the engineering and optimization of the bacterial sensor kinetics. This necessitates the development of biologically meaningful computational models for bacterial sensors comprising a variety of different molecular mechanisms, which further facilitates a systematic and quantitative evaluation of optimization strategies. In this work, we dissect the detection chain of bacterial sensors from a mathematical perspective from which we derive, supported by wet-lab data, a complete computational model for a Streptococcus mutans-based bacterial sensor as a case example. We address the engineering of bacterial sensors by investigating the impact of altered bacterial cell properties on the sensor response characteristics, specifically sensor sensitivity and response signal intensity. This is achieved through a sensitivity analysis targeting both the steady-state and transient sensor response characteristics. Alongside the demonstration of suitability of our methodological approach, our analysis shows that an increase of sensor sensitivity, through a targeted manipulation of bacterial physiology, often comes at the cost of generally diminished sensor response intensity.

Omics data analysis is crucial for studying complex diseases, but its high dimensionality and heterogeneity challenge classical statistical and machine learning methods. Graph neural networks have emerged as promising alternatives, yet the optimal strategies for their design and optimization in real-world biomedical challenges remain unclear. This study evaluates various graph representation learning models for case-control classification using high-throughput biological data from Parkinson's disease and control samples. We compare topologies derived from sample similarity networks and molecular interaction networks, including protein-protein and metabolite-metabolite interactions (PPI, MMI). Graph Convolutional Network (GCNs), Chebyshev spectral graph convolution (ChebyNet), and Graph Attention Network (GAT), are evaluated alongside advanced architectures like graph transformers, the graph U-net, and simpler models like multilayer perceptron (MLP). These models are systematically applied to transcriptomics and metabolomics data independently. Our comparative analysis highlights the benefits and limitations of various architectures in extracting patterns from omics data, paving the way for more accurate and interpretable models in biomedical research., Comment: Submitted to Machine Learning in Computational Biology 2024 as an extended abstract, 2 pages + 1 appendix

Determining gene regulatory network (GRN) structure is a central problem in biology, with a variety of inference methods available for different types of data. For a widely prevalent and challenging use case, namely single-cell gene expression data measured after intervention at multiple time points with unknown joint distributions, there is only one known specifically developed method, which does not fully utilize the rich information contained in this data type. We develop an inference method for the GRN in this case, netWork infErence by covariaNce DYnamics, dubbed WENDY. The core idea of WENDY is to model the dynamics of the covariance matrix, and solve this dynamics as an optimization problem to determine the regulatory relationships. To evaluate its effectiveness, we compare WENDY with other inference methods using synthetic data and experimental data. Our results demonstrate that WENDY performs well across different data sets.

In this paper, a set of Python methods is described that can be used to compute the frequency response of an arbitrary biochemical network given any input and output. Models can be provided in standard SBML or Antimony format. The code takes into account any conserved moieties so that this software can be used to also study signaling networks where moiety cycles are common. A utility method is also provided to make it easy to plot standard Bode plots from the generated results. The code also takes into account the possibility that the phase shift could exceed 180 degrees which can result in ugly discontinuities in the Bode plot. In the paper, some of the theory behind the method is provided as well as some commentary on the code and several illustrative examples to show the code in operation. Illustrative examples include linear reaction chains of varying lengths and the effect of negative feedback on the frequency response. Software License: MIT Open Source Availability: The code is available from https://github.com/sys-bio/frequencyResponse., Comment: Fixed some minor typos

The Hilbert number $H(n)$ is defined as the maximum number of limit cycles of a planar autonomous system of ordinary differential equations (ODEs) with right-hand sides containing polynomials of degree at most $n \in {\mathbb N}$. The dynamics of chemical reaction systems with two chemical species can be (under mass-action kinetics) described by such planar autonomous ODEs, where $n$ is equal to the maximum order of the chemical reactions in the system. Generalizations of the Hilbert number $H(n)$ to three different classes of chemical reaction networks are investigated: (i) chemical systems with reactions up to the $n$-th order; (ii) systems with up to $n$-molecular chemical reactions; and (iii) weakly reversible chemical reaction networks. In each case (i), (ii) and (iii), the question on the number of limit cycles is considered. Lower bounds on the generalized Hilbert numbers are provided for both algebraic and non-algebraic limit cycles. Furthermore, given a general algebraic curve $h(x,y)=0$ of degree $n_h \in {\mathbb N}$ and containing one or more ovals in the positive quadrant, a chemical system is constructed which has the oval(s) as its stable algebraic limit cycle(s). The ODEs describing the dynamics of the constructed chemical system contain polynomials of degree at most $n=2\,n_h+1.$ Considering $n_h \ge 4,$ the algebraic curve $h(x,y)=0$ can contain multiple closed components with the maximum number of ovals given by Harnack's curve theorem as $1+(n_h-1)(n_h-2)/2$, which is equal to 4 for $n_h=4.$ Algebraic curve $h(x,y)=0$ with $n_h=4$ and the maximum number of four ovals is used to construct a chemical system which has four stable algebraic limit cycles.

Dynamical systems with polynomials on the right-hand side can model a wide range of physical processes. A subset of such dynamical systems that can model chemical reactions under mass-action kinetics are called chemical systems. A central problem in synthetic biology is to map general polynomial dynamical systems into dynamically similar chemical ones. In this paper, we present a novel map, called the quasi-chemical map, that can systematically solve this problem. The quasi-chemical map introduces suitable state-dependent perturbations into any given polynomial dynamical system which then becomes chemical under suitably large translation of variables. We prove that this map preserves robust dynamical features, such as generic equilibria and limit cycles, as well as temporal properties, such as periods of oscillations. Furthermore, the resulting chemical systems are of only at most one degree higher than the original dynamical systems. We demonstrate the quasi-chemical map by designing relatively simple chemical systems with exotic dynamics and pre-defined bifurcation structures.

Many important phenomena in chemistry and biology are realized via dynamical features such as multi-stability, oscillations, and chaos. Construction of novel chemical systems with such finely-tuned dynamics is a challenging problem central to the growing field of synthetic biology. In this paper, we address this problem by putting forward a molecular version of a recurrent artificial neural network, which we call a recurrent neural chemical reaction network (RNCRN). We prove that the RNCRN, with sufficiently many auxiliary chemical species and suitable fast reactions, can be systematically trained to achieve any dynamics. This approximation ability is shown to hold independent of the initial conditions for the auxiliary species, making the RNCRN more experimentally feasible. To demonstrate the results, we present a number of relatively simple RNCRNs trained to display a variety of biologically-important dynamical features.

We review the construction and evolution of mathematical models of the Arabidopsis circadian clock, structuring the discussion into two distinct historical phases of modeling strategies: extension and reduction. The extension phase explores the bottom-up assembly of regulatory networks introducing as many components and interactions as possible in order to capture the oscillatory nature of the clock. The reduction phase deals with functional decomposition, distilling complex models to their essential dynamical repertoire. Current challenges in this field, including the integration of spatial considerations and environmental influences like light and temperature, are also discussed. The review emphasizes the ongoing need for models that balance molecular detail with practical simplicity.

We demonstrate that, for a chemical reaction network (CRN) engaged in energy transduction, its optimal operation from a thermodynamic efficiency standpoint is contingent upon its working conditions. Analogously to the bicycle gear system, CRNs have at their disposal several transducing mechanisms characterized by different yields. We highlight the critical role of the CRN's elementary flux modes in determining this "gearing" and their impact on maximizing energy transduction efficiency. Furthermore, we introduce an enzymatically regulated CRN, engineered to autonomously adjust its "gear", thereby optimizing its efficiency under different external conditions., Comment: 6 pages, 3 figures

This paper presents a methodology for enhancing relation extraction from biomedical texts, focusing specifically on chemical-gene interactions. Leveraging the BioBERT model and a multi-layer fully connected network architecture, our approach integrates the ChemProt and DrugProt datasets using a novel merging strategy. Through extensive experimentation, we demonstrate significant performance improvements, particularly in CPR groups shared between the datasets. The findings underscore the importance of dataset merging in augmenting sample counts and improving model accuracy. Moreover, the study highlights the potential of automated information extraction in biomedical research and clinical practice.

Phosphorylation networks, representing the mechanisms by which proteins are phosphorylated at one or multiple sites, are ubiquitous in cell signalling and display rich dynamics such as unlimited multistability. Dual-site phosphorylation networks are known to exhibit oscillations in the form of periodic trajectories, when phosphorylation and dephosphorylation occurs as a mixed mechanism: phosphorylation of the two sites requires one encounter of the kinase, while dephosphorylation of the two sites requires two encounters with the phosphatase. A still open question is whether a mechanism requiring two encounters for both phosphorylation and dephosphorylation also admits oscillations. In this work we provide evidence in favor of the absence of oscillations of this network by precluding Hopf bifurcations in any reduced network comprising three out of its four intermediate protein complexes. Our argument relies on a novel network reduction step that preserves the absence of Hopf bifurcations, and on a detailed analysis of the semi-algebraic conditions precluding Hopf bifurcations obtained from Hurwitz determinants of the characteristic polynomial of the Jacobian of the system. We conjecture that the removal of certain reverse reactions appearing in Michaelis-Menten-type mechanisms does not have an impact on the presence or absence of Hopf bifurcations. We prove an implication of the conjecture under certain favorable scenarios and support the conjecture with additional example-based evidence.

Biological cells rely on precise spatiotemporal coordination of biochemical reactions to control their many functions. Such cell signaling networks have been a common focus for mathematical models, but they remain challenging to simulate, particularly in realistic cell geometries. Herein, we present our software, Spatial Modeling Algorithms for Reactions and Transport (SMART), a package that takes in high-level user specifications about cell signaling networks and molecular transport, and then assembles and solves the associated mathematical and computational systems. SMART uses state-of-the-art finite element analysis, via the FEniCS Project software, to efficiently and accurately resolve cell signaling events over discretized cellular and subcellular geometries. We demonstrate its application to several different biological systems, including YAP/TAZ mechanotransduction, calcium signaling in neurons and cardiomyocytes, and ATP generation in mitochondria. Throughout, we utilize experimentally-derived realistic cellular geometries represented by well-conditioned tetrahedral meshes. These scenarios demonstrate the applicability, flexibility, accuracy and efficiency of SMART across a range of temporal and spatial scales.

Antimony is a high-level, human-readable text-based language designed for defining and sharing models in the systems biology community. It enables scientists to describe biochemical networks and systems using a simple and intuitive syntax. It allows users to easily create, modify, and distribute reproducible computational models. By allowing the concise representation of complex biological processes, Antimony enhances collaborative efforts, improves reproducibility, and accelerates the iterative development of models in systems biology. This paper provides an update to the Antimony language since it was introduced in 2009. In particular, we highlight new annotation features, support for flux balance analysis, a new rateOf method, support for probability distributions and uncertainty, named stochiometries, and algebraic rules. Antimony is also now distributed as a C/C++ library, together with python and Julia bindings, as well as a JavaScript version for use within a web browser. Availability: https://github.com/sys-bio/antimony.

Drug response prediction (DRP) is a crucial phase in drug discovery, and the most important metric for its evaluation is the IC50 score. DRP results are heavily dependent on the quality of the generated molecules. Existing molecule generation methods typically employ classifier-based guidance, enabling sampling within the IC50 classification range. However, these methods fail to ensure the sampling space range's effectiveness, generating numerous ineffective molecules. Through experimental and theoretical study, we hypothesize that conditional generation based on the target IC50 score can obtain a more effective sampling space. As a result, we introduce regressor-free guidance molecule generation to ensure sampling within a more effective space and support DRP. Regressor-free guidance combines a diffusion model's score estimation with a regression controller model's gradient based on number labels. To effectively map regression labels between drugs and cell lines, we design a common-sense numerical knowledge graph that constrains the order of text representations. Experimental results on the real-world dataset for the DRP task demonstrate our method's effectiveness in drug discovery. The code is available at:https://anonymous.4open.science/r/RMCD-DBD1., Comment: 22 pages, 7 figures, 9 tables

The upcoming phase of space exploration not only includes trips to Mars and beyond, but also holds great promise for human progress. However, the harm caused by cosmic radiation, consisting of Galactic Cosmic Rays and Solar Particle Events, is an important safety concern for astronauts and other living things that will accompany them. Research exploring the biological effects of cosmic radiation includes experiments conducted in space itself and in simulated space environments on Earth. Notably, NASA's Space Radiation Laboratory has taken significant steps forward in simulating cosmic radiation by using particle accelerators and is currently pioneering the progress in this field. Curiously, much of the research emphasis thus far has been on understanding how cosmic radiation impacts living organisms, instead of finding ways to help them resist the radiation. In this paper, we briefly talk about current research on the biological effects of cosmic radiation and propose possible protective measures through biological interventions. In our opinion, biological response pathways responsible for coping with stressors on Earth can provide effective solutions for protection against the stress caused by cosmic radiation. We also recommend evaluating the effectiveness of these pathways through experiments using particle accelerators to simulate the effects of cosmic radiation.

Extant life contains numerous non-equilibrium mechanisms to create order not achievable at equilibrium; it is generally assumed that these mechanisms evolved because the resulting order was sufficiently beneficial to overcome associated costs of time and energy. Here, we identify a broad range of conditions under which non-equilibrium order-creating mechanisms will evolve as an inevitable consequence of self-replication, even if the order is not directly functional. We show that models of polymerases, when expanded to include known stalling effects, can evolve kinetic proofreading through selection for fast replication alone, consistent with data from recent mutational screens. Similarly, replication contingent on fast self-assembly can select for non-equilibrium instabilities and result in more ordered structures without any direct selection for order. We abstract these results into a framework that predicts that self-replication intrinsically amplifies dissipative order-enhancing mechanisms if the distribution of replication times is wide enough. Our work suggests the intriguing possibility that non-equilibrium order can arise more easily than assumed, even before that order is directly functional, with consequences impacting such diverse phenomena as the evolution of mutation rates, kinetic traps in self-assembly, and the origin of life.

Biomarkers enable objective monitoring of a given cell or state in a biological system and are widely used in research, biomanufacturing, and clinical practice. However, identifying appropriate biomarkers that are both robustly measurable and capture a state accurately remains challenging. We present a framework for biomarker identification based upon observability guided sensor selection. Our methods, Dynamic Sensor Selection (DSS) and Structure-Guided Sensor Selection (SGSS), utilize temporal models and experimental data, offering a template for applying observability theory to data from biological systems. Unlike conventional methods that assume well-known, fixed dynamics, DSS adaptively select biomarkers or sensors that maximize observability while accounting for the time-varying nature of biological systems. Additionally, SGSS incorporates structural information and diverse data to identify sensors which are resilient against inaccuracies in our model of the underlying system. We validate our approaches by performing estimation on high dimensional systems derived from temporal gene expression data from partial observations. Our algorithms reliably identify known biomarkers and uncover new ones within our datasets. Additionally, integrating chromosome conformation and gene expression data addresses noise and uncertainty, enhancing the reliability of our biomarker selection approach for the genome.

Understanding the insulin signaling cascade provides insights on the underlying mechanisms of biological phenomena such as insulin resistance, diabetes, Alzheimer's disease, and cancer. For this reason, previous studies utilized chemical reaction network theory to perform comparative analyses of reaction networks of insulin signaling in healthy (INSMS: INSulin Metabolic Signaling) and diabetic cells (INRES: INsulin RESistance). This study extends these analyses using various methods which give further insights regarding insulin signaling. Using embedded networks, we discuss evidence of the presence of a structural "bifurcation" in the signaling process between INSMS and INRES. Concordance profiles of INSMS and INRES show that both have a high propensity to remain monostationary. Moreover, the concordance properties allow us to present heuristic evidence that INRES has a higher level of stability beyond its monostationarity. Finally, we discuss a new way of analyzing reaction networks through network translation. This method gives rise to three new insights: (i) each stoichiometric class of INSMS and INRES contains a unique positive equilibrium; (ii) any positive equilibrium of INSMS is exponentially stable and is a global attractor in its stoichiometric class; and (iii) any positive equilibrium of INRES is locally asymptotically stable. These results open up opportunities for collaboration with experimental biologists to understand insulin signaling better., Comment: 18 pages, 0 figure

Proteins congregate into complexes to perform fundamental cellular functions. Phenotypic outcomes, in health and disease, are often mechanistically driven by the remodeling of protein complexes by protein coding mutations or cellular signaling changes in response to molecular cues. Here, we present an affinity purification mass spectrometry (APMS) proteomics protocol to quantify and visualize global changes in protein protein interaction (PPI) networks between pairwise conditions. We describe steps for expressing affinity tagged bait proteins in mammalian cells, identifying purified protein complexes, quantifying differential PPIs, and visualizing differential PPI networks. Specifically, this protocol details steps for designing affinity tagged bait gene constructs, transfection, affinity purification, mass spectrometry sample preparation, data acquisition, database search, data quality control, PPI confidence scoring, cross run normalization, statistical data analysis, and differential PPI visualization. Our protocol discusses caveats and limitations with applicability across cell types and biological areas., Comment: 29 pages, 3 figures

Drug development is a lengthy process with a high failure rate. Increasingly, machine learning is utilized to facilitate the drug development processes. These models aim to enhance our understanding of drug characteristics, including their activity in biological contexts. However, a major challenge in drug response (DR) prediction is model interpretability as it aids in the validation of findings. This is important in biomedicine, where models need to be understandable in comparison with established knowledge of drug interactions with proteins. drGAT, a graph deep learning model, leverages a heterogeneous graph composed of relationships between proteins, cell lines, and drugs. drGAT is designed with two objectives: DR prediction as a binary sensitivity prediction and elucidation of drug mechanism from attention coefficients. drGAT has demonstrated superior performance over existing models, achieving 78\% accuracy (and precision), and 76\% F1 score for 269 DNA-damaging compounds of the NCI60 drug response dataset. To assess the model's interpretability, we conducted a review of drug-gene co-occurrences in Pubmed abstracts in comparison to the top 5 genes with the highest attention coefficients for each drug. We also examined whether known relationships were retained in the model by inspecting the neighborhoods of topoisomerase-related drugs. For example, our model retained TOP1 as a highly weighted predictive feature for irinotecan and topotecan, in addition to other genes that could potentially be regulators of the drugs. Our method can be used to accurately predict sensitivity to drugs and may be useful in the identification of biomarkers relating to the treatment of cancer patients.

The biomolecular networks underpinning cell function exhibit canalization, or the buffering of fluctuations required to function in a noisy environment. One understudied putative mechanism for canalization is the functional equivalence of a biomolecular entity's regulators (e.g., among the transcription factors for a gene). In these discrete dynamical systems, activation and inhibition of biomolecular entities (e.g., transcription of genes) are modeled as the activity of coupled 2-state automata, and thus the equivalence of regulators can be studied using the theory of symmetry in discrete functions. To this end, we present a new exact algorithm for finding maximal symmetry groups among the inputs to discrete functions. We implement this algorithm in Rust as a Python package, schematodes. We include schematodes in the new CANA v1.0.0 release, an open source Python library for analyzing canalization in Boolean networks, which we also present here. We compare our exact method implemented in schematodes to the previously published inexact method used in earlier releases of CANA and find that schematodes significantly outperforms the prior method both in speed and accuracy. We also apply CANA v1.0.0 to study the symmetry properties of regulatory function from an ensemble of experimentally-supported Boolean networks from the Cell Collective. Using CANA v1.0.0, we find that the distribution of a previously reported symmetry parameter, $k_s/k$, is statistically significantly different in the Cell Collective than in random automata with the same in-degree and activation bias (Kolmogorov-Smirnov test, $p<0.001$). In particular, its spread is much wider than in our null model (IQR 0.31 vs IQR 0.20 with equal medians), demonstrating that the Cell Collective is enriched in functions with extreme symmetry or asymmetry., Comment: 5 pages 1 figure (two images in figure)