Distilling dynamical knowledge from stochastic reaction networks.

Stochastic reaction networks are widely used in the modeling of stochastic systems across diverse domains such as biology, chemistry, physics, and ecology. However, the comprehension of the dynamic behaviors inherent in stochastic reaction networks is a formidable undertaking, primarily due to the exponential growth in the number of possible states or trajectories as the state space dimension increases. In this study, we introduce a knowledge distillation method based on reinforcement learning principles, aimed at compressing the dynamical knowledge encoded in stochastic reaction networks into a singular neural network construct. The trained neural network possesses the capability to accurately predict the state conditional joint probability distribution that corresponds to the given query contexts, when prompted with rate parameters, initial conditions, and time values. This obviates the need to track the dynamical process, enabling the direct estimation of normalized state and trajectory probabilities, without necessitating the integration over the complete state space. By applying our method to representative examples, we have observed a high degree of accuracy in both multimodal and high-dimensional systems. Additionally, the trained neural network can serve as a foundational model for developing efficient algorithms for parameter inference and trajectory ensemble generation. These results collectively underscore the efficacy of our approach as a universal means of distilling knowledge from stochastic reaction networks. Importantly, our methodology also spotlights the potential utility in harnessing a singular, pretrained, large-scale model to encapsulate the solution space underpinning a wide spectrum of stochastic dynamical systems.

A mixed-effects stochastic model reveals clonal dominance in gene therapy safety studies.

BACKGROUND: Mathematical models of haematopoiesis can provide insights on abnormal cell expansions (clonal dominance), and in turn can guide safety monitoring in gene therapy clinical applications. Clonal tracking is a recent high-throughput technology that can be used to quantify cells arising from a single haematopoietic stem cell ancestor after a gene therapy treatment. Thus, clonal tracking data can be used to calibrate the stochastic differential equations describing clonal population dynamics and hierarchical relationships in vivo. RESULTS: In this work we propose a random-effects stochastic framework that allows to investigate the presence of events of clonal dominance from high-dimensional clonal tracking data. Our framework is based on the combination between stochastic reaction networks and mixed-effects generalized linear models. Starting from the Kramers-Moyal approximated Master equation, the dynamics of cells duplication, death and differentiation at clonal level, can be described by a local linear approximation. The parameters of this formulation, which are inferred using a maximum likelihood approach, are assumed to be shared across the clones and are not sufficient to describe situation in which clones exhibit heterogeneity in their fitness that can lead to clonal dominance. In order to overcome this limitation, we extend the base model by introducing random-effects for the clonal parameters. This extended formulation is calibrated to the clonal data using a tailor-made expectation-maximization algorithm. We also provide the companion  package RestoreNet, publicly available for download at https://cran.r-project.org/package=RestoreNet . CONCLUSIONS: Simulation studies show that our proposed method outperforms the state-of-the-art. The application of our method in two in-vivo studies unveils the dynamics of clonal dominance. Our tool can provide statistical support to biologists in gene therapy safety analyses.

Stochastic approximations of higher-molecular by bi-molecular reactions.

Reactions involving three or more reactants, called higher-molecular reactions, play an important role in mathematical modelling in systems and synthetic biology. In particular, such reactions underpin a variety of important bio-dynamical phenomena, such as multi-stability/multi-modality, oscillations, bifurcations, and noise-induced effects. However, as opposed to reactions involving at most two reactants, called bi-molecular reactions, higher-molecular reactions are biochemically improbable. To bridge the gap, in this paper we put forward an algorithm for systematically approximating arbitrary higher-molecular reactions with bi-molecular ones, while preserving the underlying stochastic dynamics. Properties of the algorithm and convergence are established via singular perturbation theory. The algorithm is applied to a variety of higher-molecular biochemical networks, and is shown to play an important role in synthetic biology.

Stationary distributions via decomposition of stochastic reaction networks.

We examine reaction networks (CRNs) through their associated continuous-time Markov processes. Studying the dynamics of such networks is in general hard, both analytically and by simulation. In particular, stationary distributions of stochastic reaction networks are only known in some cases. We analyze class properties of the underlying continuous-time Markov chain of CRNs under the operation of join and examine conditions such that the form of the stationary distributions of a CRN is derived from the parts of the decomposed CRNs. The conditions can be easily checked in examples and allow recursive application. The theory developed enables sequential decomposition of the Markov processes and calculations of stationary distributions. Since the class of processes expressible through such networks is big and only few assumptions are made, the principle also applies to other stochastic models. We give examples of interest from CRN theory to highlight the decomposition.

Sensitivity Analysis for Multiscale Stochastic Reaction Networks Using Hybrid Approximations.

We consider the problem of estimating parameter sensitivities for stochastic models of multiscale reaction networks. These sensitivity values are important for model analysis, and the methods that currently exist for sensitivity estimation mostly rely on simulations of the stochastic dynamics. This is problematic because these simulations become computationally infeasible for multiscale networks due to reactions firing at several different timescales. However it is often possible to exploit the multiscale property to derive a "model reduction" and approximate the dynamics as a Piecewise deterministic Markov process, which is a hybrid process consisting of both discrete and continuous components. The aim of this paper is to show that such PDMP approximations can be used to accurately and efficiently estimate the parameter sensitivity for the original multiscale stochastic model. We prove the convergence of the original sensitivity to the corresponding PDMP sensitivity, in the limit where the PDMP approximation becomes exact. Moreover, we establish a representation of the PDMP parameter sensitivity that separates the contributions of discrete and continuous components in the dynamics and allows one to efficiently estimate both contributions.

Non-explosivity of Stochastically Modeled Reaction Networks that are Complex Balanced.

We consider stochastically modeled reaction networks and prove that if a constant solution to the Kolmogorov forward equation decays fast enough relatively to the transition rates, then the model is non-explosive. In particular, complex-balanced reaction networks are non-explosive.

A stochastic transcriptional switch model for single cell imaging data.

Gene expression is made up of inherently stochastic processes within single cells and can be modeled through stochastic reaction networks (SRNs). In particular, SRNs capture the features of intrinsic variability arising from intracellular biochemical processes. We extend current models for gene expression to allow the transcriptional process within an SRN to follow a random step or switch function which may be estimated using reversible jump Markov chain Monte Carlo (MCMC). This stochastic switch model provides a generic framework to capture many different dynamic features observed in single cell gene expression. Inference for such SRNs is challenging due to the intractability of the transition densities. We derive a model-specific birth-death approximation and study its use for inference in comparison with the linear noise approximation where both approximations are considered within the unifying framework of state-space models. The methodology is applied to synthetic as well as experimental single cell imaging data measuring expression of the human prolactin gene in pituitary cells.

Fast reactions with non-interacting species in stochastic reaction networks.

We consider stochastic reaction networks modeled by continuous-time Markov chains. Such reaction networks often contain many reactions, potentially occurring at different time scales, and have unknown parameters (kinetic rates, total amounts). This makes their analysis complex. We examine stochastic reaction networks with non-interacting species that often appear in examples of interest (e.g. in the two-substrate Michaelis Menten mechanism). Non-interacting species typically appear as intermediate (or transient) chemical complexes that are depleted at a fast rate. We embed the Markov process of the reaction network into a one-parameter family under a two time-scale approach, such that molecules of non-interacting species are degraded fast. We derive simplified reaction networks where the non-interacting species are eliminated and that approximate the scaled Markov process in the limit as the parameter becomes small. Then, we derive sufficient conditions for such reductions based on the reaction network structure for both homogeneous and time-varying stochastic settings, and study examples and properties of the reduction.

Antithetic proportional-integral feedback for reduced variance and improved control performance of stochastic reaction networks.

The ability of a cell to regulate and adapt its internal state in response to unpredictable environmental changes is called homeostasis and this ability is crucial for the cell's survival and proper functioning. Understanding how cells can achieve homeostasis, despite the intrinsic noise or randomness in their dynamics, is fundamentally important for both systems and synthetic biology. In this context, a significant development is the proposed antithetic integral feedback (AIF) motif, which is found in natural systems, and is known to ensure robust perfect adaptation for the mean dynamics of a given molecular species involved in a complex stochastic biomolecular reaction network. From the standpoint of applications, one drawback of this motif is that it often leads to an increased cell-to-cell heterogeneity or variance when compared to a constitutive (i.e. open-loop) control strategy. Our goal in this paper is to show that this performance deterioration can be countered by combining the AIF motif and a negative feedback strategy. Using a tailored moment closure method, we derive approximate expressions for the stationary variance for the controlled network that demonstrate that increasing the strength of the negative feedback can indeed decrease the variance, sometimes even below its constitutive level. Numerical results verify the accuracy of these results and we illustrate them by considering three biomolecular networks with two types of negative feedback strategies. Our computational analysis indicates that there is a trade-off between the speed of the settling-time of the mean trajectories and the stationary variance of the controlled species; i.e. smaller variance is associated with larger settling-time.

Adaptive moment closure for parameter inference of biochemical reaction networks.

Continuous-time Markov chain (CTMC) models have become a central tool for understanding the dynamics of complex reaction networks and the importance of stochasticity in the underlying biochemical processes. When such models are employed to answer questions in applications, in order to ensure that the model provides a sufficiently accurate representation of the real system, it is of vital importance that the model parameters are inferred from real measured data. This, however, is often a formidable task and all of the existing methods fail in one case or the other, usually because the underlying CTMC model is high-dimensional and computationally difficult to analyze. The parameter inference methods that tend to scale best in the dimension of the CTMC are based on so-called moment closure approximations. However, there exists a large number of different moment closure approximations and it is typically hard to say a priori which of the approximations is the most suitable for the inference procedure. Here, we propose a moment-based parameter inference method that automatically chooses the most appropriate moment closure method. Accordingly, contrary to existing methods, the user is not required to be experienced in moment closure techniques. In addition to that, our method adaptively changes the approximation during the parameter inference to ensure that always the best approximation is used, even in cases where different approximations are best in different regions of the parameter space.

Solution of Chemical Master Equations for Nonlinear Stochastic Reaction Networks.

Stochasticity in the dynamics of small reacting systems requires discrete-probabilistic models of reaction kinetics instead of traditional continuous-deterministic ones. The master probability equation is a complete model of randomly evolving molecular populations. Because of its ambitious character, the master equation remained unsolved for all but the simplest of molecular interaction networks. With the first solution of chemical master equations, a wide range of experimental observations of small-system interactions may be mathematically conceptualized.

An efficient and unbiased method for sensitivity analysis of stochastic reaction networks.

We consider the problem of estimating parameter sensitivity for Markovian models of reaction networks. Sensitivity values measure the responsiveness of an output with respect to the model parameters. They help in analysing the network, understanding its robustness properties and identifying the important reactions for a specific output. Sensitivity values are commonly estimated using methods that perform finite-difference computations along with Monte Carlo simulations of the reaction dynamics. These methods are computationally efficient and easy to implement, but they produce a biased estimate which can be unreliable for certain applications. Moreover, the size of the bias is generally unknown and hence the accuracy of these methods cannot be easily determined. There also exist unbiased schemes for sensitivity estimation but these schemes can be computationally infeasible, even for very simple networks. Our goal in this paper is to present a new method for sensitivity estimation, which combines the computational efficiency of finite-difference methods with the accuracy of unbiased schemes. Our method is easy to implement and it relies on an exact representation of parameter sensitivity that we recently proved elsewhere. Through examples, we demonstrate that the proposed method can outperform the existing methods, both biased and unbiased, in many situations.