article: The dynamics of stochastic reaction networks within cells are inevitably modulated by factors considered extrinsic to the network such as, for instance, the fluctuations in ribosome copy numbers for a gene regulatory network. While several recent studies demonstrate the importance of accounting for such extrinsic components, the resulting models are typically hard to analyze. In this work we develop a general mathematical framework that allows to uncouple the network from its dynamic environment by incorporating only the environment' s effect onto the network into a new model. More technically, we show how such fluctuating extrinsic components (e. g. , chemical species) can be marginalized in order to obtain this decoupled model. We derive its corresponding process- and master equations and show how stochastic simulations can be performed. Using several case studies, we demonstrate the significance of the approach.
  Biochemical systems involving low-copy molecules demand for mathematical models that account for the intrinsic stochasticity [1]. In recent years, however, realization has grown that intrinsic noise alone cannot account for the observed substantial phenotypic variability among isogenic cells. That is, fluctuations in the intracellular environment, commonly termed extrinsic noise, represent an additional source of variability [2], [3], [4]. Several recent studies focus on separating intrinsic and extrinsic fluctuations through dual-reporter measurements [5], [6]. Other approaches model extrinsic noise through certain parameters (i. e. , the translation rate) of a kinetic model which is calibrated subsequently using flow-cytometry [7], [8] or time-lapse microscopy data [9], [10]. All of those approaches have in common that they consider the biochemical process under study – i. e. , the expression of a gene – as a small subpart that is embedded into a larger dynamical system. Accordingly, they rely on augmenting the original kinetic model by certain environmental components which are assumed to be fixed but random [11], [10], [8], [9] or fluctuating over time [5], [12], [13]. In fact, such models agree very well with the variability that is observed experimentally, but on their downside, suffer from the increased dimensionality - somehow defeating the original purpose of tractable dedicated models. A natural question arising in that context is whether we can find a proper dynamical description of just the system of interest as if it was still embedded into its stochastically modulating environment. In other words, we aim to find a “self-contained” stochastic model that summarizes all system behaviors attainable under all possible realizations of the extrinsic fluctuations. Such models could then be used to perform an uncoupled analysis of a reaction network subject to extrinsic noise. The mathematical correct answer that we provide in this work is the marginalization of the system dynamics with respect to those extrinsic fluctuations. Interestingly it turns out that the resulting model exploits its own stochasticity to emulate the effect of extrinsic noise, leading to a self-exciting process. A simple instance of such self-excitation is the Polya urn scheme: at each draw from the Polya urn with balls of two colors the drawn ball and a fixed number of new balls of the same color as the draw are placed in the urn. This scheme is known to be equivalent to Bernoulli trials marginalized over random (and here correspondingly extrinsic) success rates [14]. Intuitively, due to its self-excitation, the number of draws of the same color for a Polya urn over repeated draws displays a much richer and dispersed dynamics than the number of successful draws in a Bernoulli trial with fixed success rate. For the purpose of inference we recently proposed a first attempt of such marginalization for the special case of fixed but random environmental conditions [9]. In this work we develop a general mathematical framework from which the uncoupled dynamics can be constructed in a principled manner, regardless whether the environment is constant or dynamically changing.
  We describe the time-evolution of a stochastic reaction network by a continous-time Markov chain (CTMC) with chemical species and reaction channels. The system state at time is denoted and we write its random path on time intervals as. Throughout we will follow the usual convention to refer to upper-case and lower-case versions of a symbol as a random variable and its realization, respectively. Furthermore, we assume that depends on another multivariate Markov process through its hazard functions in the form (1) with some positive function and a polynomial determined by the law of mass-action, for instance. For reactions independent of, we thus have. Typically, is another jump or diffusion process corresponding to a set of modulating environmental species or conditions that are considered extrinsic to the system of interest, whereas the species in represent the actual system of interest. For example, could be the fluctuating ribosome copy numbers affecting the kinetics of a gene regulatory network represented by. Although a more general treatment is possible, we assume a feed-forward structure between and, which means that modulates but not vice-versa. Consequently, the dynamics of the joint system can be described by a Markov process together with a conditional Markov chain.
  There is increasing evidence that models of biochemical networks need to account for both intrinsic and extrinsic noise caused by variations in the intracellular environment. In recent studies, this is done by extending a model' s state space by certain environmental species, whose dynamics are described along with the actual system of interest. In particular, the resulting system dynamics are described and studied conditional on a particular history of the environment and thus, do not provide a coherent description of a dynamical system subject to extrinsic noise. In this work, we derived and analyzed a novel process framework, which is able to describe just the system of interest as if it was still embedded into its environment. In that sense, it permits a mathematically exact way to analyze small parts of networks in an uncoupled fashion. Several recent studies rely on the extreme assumptions that the environmental fluctuations are either infinitely fast or slow. While both strategies may in fact lead to strongly simplified and tractable models, they are characterized by significant approximation errors when considering intermediate environmental timescales (see e. g. , Fig. 5b). The approach proposed here allows to uncouple a reaction network from its surrounding environment regardless of the latter' s timescale. In that sense, the approach is fully general although practical implementations may rely on efficient but approximate solutions of the discussed filtering problem. In the context of Monte Carlo simulation the decoupled process can yield a significant reduction in computational effort when compared to standard SSA – especially if the environmental network is costly to simulate. This highlights the role of the provided framework as a general tool to split stochastic biochemical networks into individual parts that are easier to simulate. We believe that it will aid in turning stochastic modeling and simulation techniques more large-scale and more faithful to in vivo conditions, where significant environmental fluctuations are present. Moreover, the framework can be used in the model-based design of novel circuit motifs in synthetic biology and is related to the notion of retroactivity [35]. We further demonstrate that the uncoupled dynamics provide a novel analytical tool to study how environmental stochasticity is propagated along coupled reaction networks. For instance, we have shown that the total environmental noise splits up into two terms: one corresponding to the noise that is suppressed and a second term that quantifies the effective noise that is sensed by the target network. In [36] the authors derive a lower bound on a network' s ability to suppress fluctuations and show its immediate relation to the uncertainty at which those fluctuations can be estimated – similar to what we defined as effective noise. The methods proposed here allow to not only bound, but fully determine both the suppressed and effective environmental noise. Our results further indicate that two environments with very different timescales may impact a network in a similar way. For instance a fixed but random environment may yield the same effective noise as a fluctuating environment with larger variance. Along those lines, we derived a simple but widely applicable approximation of the transient probability distribution for birth-death processes subject to environmental noise. It is based on the idea to approximate a fluctuating environmental process by a simple random variable that impacts the birth death process in an equivalent way. In order to solve for the transient probability distribution we derived a generalized master equation for this non-Markovian process. The results were derived under the assumption of a one-dimensional environmental process but its extension to multivariate scenarios shall be examined at a later stage. The presented framework provides a novel and orthogonal notion of environmentally perturbed stochastic networks from which several interesting avenues could be pursued. However, the approach in its current form has a few practical limitations that should be addressed in future work. Above all, if the uncoupled dynamics cannot be computed analytically, its accuracy – and hence, reliability will depend on finding suitable approximate schemes. For that sake, we believe that one will profit from the approach' s tight connection to stochastic filtering – a field that has been thoroughly explored over the past decades. Therefore, the marginal hazard functions may be approximated by exploiting certain features of the associated filtering problems such as their asymptotic theories and so forth. Those could for instance be used to derive more general moment-closure schemes or to construct concise and more efficient state-space truncations [37] when integrating the respective filtering distribution. The main bottleneck of the MSA algorithm is the simulation of the non-exponential waiting-times, which relies on the integration of an ODE. Therefore, the use of optimized and dedicated ODE solvers could further boost the algorithm' s performance.
  All models and simulations were implemented in MATLAB (The MathWorks, Natick, MA). Source codes for reproducing the numerical results have been attached as S1 Database.
summary: Cellular reaction networks show a substantial degree of complexity, and many of the regulatory mechanisms are yet to be discovered. Therefore, building detailed mathematical models of complete networks is not only unrewarding but most often impossible. Instead, focus is put on small subnetworks such as the transcriptional circuitry associated with a few genes. Practically, those subnetworks are still embedded in the cellular environment and are consequently modulated by it. The mathematical framework that we present in this work allows to build models of such subnetworks only, while still incorporating the impact of the environment. Such models are instrumental not only for advancing our biological understanding through them, but also for the reliable forward engineering of new subnetworks in synthetic biology.
keywords: biochemistry
  biochemical simulations
  computer and information sciences
  genetic networks
  mathematics
  markov models
  gene regulatory networks
  genome analysis
  network analysis
  biology and life sciences
  physical sciences
  computational biology
  stochastic processes
  probability theory
title: Uncoupled Analysis of Stochastic Reaction Networks in Fluctuating Environments