Computation of stationary distributions in stochastic models of cellular processes with molecular memory

Abstract Modeling stochastic dynamics of intracellular processes has long rested on Markovian (i.e., memoryless) hypothesis. However, many of these processes are non-Markovian (i.e., memorial) due to, e.g., small reaction steps involved in synthesis or degradation of a macroscopic molecule. When interrogating aspects of a cellular network by experimental measurements (e.g., by singlemolecule and single-cell measurement technologies) of network components, a key need is to develop efficient approaches to simulate and compute joint distributions of these components. To cope with this computational challenge, we develop two efficient algorithms: stationary generalized Gillespie algorithm and stationary generalized finite state projection, both being established based on a stationary generalized chemical master equation. We show how these algorithms can be combined in a streamlined procedure for evaluation of non-Markovian effects in a general cellular network. Stationary distributions are evaluated in two models of constitutive and bursty gene expressions as well as a model of genetic toggle switch, each considering molecular memory. Our approach significantly expands the capability of stochastic simulation to investigate gene regulatory network dynamics, which has the potential to advance both understanding of molecular systems biology and design of synthetic circuits. Author summary Cellular systems are driven by interactions between subsystems via time-stamped discrete events, involving numerous components and reaction steps and spanning several time scales. Such biochemical reactions are subject to inherent noise due to the small numbers of molecules. Also, they could involve several small steps, creating a memory between individual events. Delineating these molecular stochasticity and memory of biomolecular networks are continuing challenges for molecular systems biology. We present a novel approach to compute the probability distribution in stochastic models of cellular processes with molecular memory based on stationary generalized chemical master equation. We map a stochastic system with memory onto a Markovian model with effective reaction propensity functions. This formulation enables us to efficiently develop algorithms under the Markovian framework, and thus systematically analyze how molecular memories regulate stochastic behaviors of biomolecular networks. Here we propose two representative algorithms: stationary generalized Gillespie algorithm and stationary generalized finite state projection algorithm. The former generate realizations with Monte Carlo simulation, but the later compute approximations of the probability distribution by solving a truncated version of stochastic process. Our approach is demonstrated by applying it to three different examples from systems biology: generalized birth-death process, a stochastic toggle switch model, and a 3-stage gene expression model.


Introduction
Cellular processes are biochemical and hence noisy: in each cell, the numbers or concentrations of species molecules are subject to random fluctuations due to the small numbers of these molecules or environmental perturbations or both. This stochasticity, which has been verified by singlemolecule experiments and can have a significant impact on cellular functioning, has been intensively studied in the past few decades [1,2]. For a biochemical reaction system, if the stochastic motion of the reactants is assumed to be uninfluenced by previous states and only by the current state [3,4], then the stochastic effect can be well revealed by analyzing the chemical master equation (CME) [3][4][5]. This Markovian assumption has also led to important successes in the description of many stochastic processes on networks [3][4][5][6][7][8].
However, as a general rule, the dynamics of a given reactive species results from its interactions with the environment and cannot be described as a Markov process. Indeed, although the evolution of the set of all microscopic degrees of freedom of a biochemical system is Markovian, the dynamics restricted to the species only is not [9,10] (Fig.1 (A)). This is typically the case for a promoter cycle, whose non-Markovian dynamics result from the fact that the promoter gradually becomes mature 3 through a series of inactive states (with several hidden Poisson steps) before being activated, leading to narrowly distributed gestation periods between transcription windows [11]. More generally, the expressions of most genes are controlled in a complex manner involving several repressors, transcription factors and mediators, as well as chromatin remodeling or changes in supercoiling.
These unspecified regulations or other processes can be modelled by non-exponential waiting time distributions [11][12][13][14][15]. Other experimentally observed examples of non-Markovian dynamics include the human activity patterns [16], natural phenomena [17], biological processes [18], biochemical reactions [19,20], etc. Examples listed in Table 1 indicate that waiting-time distributions between biochemical events are not limited to exponential distributions but may be diverse. Indeed, the time scales of ultrafast biochemical processes are of the same order as the time scale of the bath they interact with, leading to failure of the Markov approximation. An important step towards a general understanding of non-Markovian dynamics was recently made based on the continuous time random walk (CTRW) theory. Namely, a generalized chemical master equation (gCME) was derived, which can account for non-exponential inter-reaction times and the resulting non-Markovian character of reaction dynamics in time [26,27]. However, this equation is rather formal, and the "memory functions" involved are implicitly expressed by waiting time distributions. This leads to notorious difficulties in analysis and simulation, greatly limiting applications of the gCME. An alternative method is the kinetic Monte Carlo approach, which instead seeks either exact or approximate realizations of the underlying non-Markov process. Earlier, along the line of CTRW, Gillespie developed a stochastic simulation algorithm of random walk with residence time depending on transition probability rates [28], which is actually a prototype of the famous Gillespie algorithm (GA) [29]. Based on the renewal theory, Boguñá and colleagues 4 extended the GA to general renewal processes [30]. Subsequently, Vestergaard and Génois developed a temporal GA, which can simulate non-Markovian dynamics [31]. Recently, Masuda and Rocha proposed a Laplace GA, which can be applied to renewal processes where survival functions are monotonic [32]. In a word, these numerical methods generate realizations enough to obtain statistics for pivotal events. Unfortunately, in many cases, biologically important events occur rarely when reaction waiting times span multiple time scales or/and follow non-exponential distributions, implying that the previous methods are unfeasible in this case. It is an urgent need but is still a significant challenge to develop effective methods to simulate and compute the stochastic dynamics of non-Markovian reaction processes.
As a fact of matter, the emergence of measurement techniques such as single-cell fluorescent microscopy [33], flow and mass cytometry [34], single-cell qPCR [35] and single-cell RNA-seq [36] currently allow for access to increasingly rich data on approximately steady-state distributions of interested molecular species and waiting time distributions of key reaction events. Indeed, steadystate distribution is an important quantity needed to characterize the stationary behavior of many biochemical network systems. To assess how experimental data can be informative, it is also needed to calculate or simulate aspects of steady-state distribution [37][38][39][40][41].  Here we propose two efficient stochastic simulation algorithms: stationary generalized GA (sgGA) that generates realizations with Monte Carlo simulation, and stationary generalized finite state projection (sgFSP) that computes approximate probability distributions. These algorithms are established based on a stationary gCME (sgCME) for a general biomolecular network with arbitrary waiting time distributions. While a key point of the sgCME is that it converts a computationally challenging non-Markovian problem into a computable Markovian one ( Fig.1 (B)) by introducing an effective transition rate for every reaction in the network [27], our algorithms are designed to sufficiently exploit this characteristic and hence can be used not only in simulation and analysis of non-Markovian dynamics but also in discovery of biological knowledge underlying non-Markovian effects in complex cellular network systems.
The paper is organized as follows. First, we introduce a general cellular network model with any waiting-time distributions between biochemical reactions. Then, we model this network with a gCME expressed in Laplace transforms and derive a sgCME where effective transition rates decoding non-Markovian effects are explicitly given. Next, we propose two efficient algorithms: sgGA and sgFSP. As an illustration, we have applied these methods to the analysis of a wellcharacterized birth-death process and investigate the memory effect on stochastic gene expression.
Furthermore, we also demonstrated the applicability of these methods to two more complex biochemical network: a stochastic toggle switch model, and a 3-stage stochastic gene expression model. Finally, we revisit derivation of the most common waiting time distributions in reaction systems and give their physical explanations, with the aim to help the reader's understanding.

A general theory for reaction networks
Consider a reaction network that models a general cellular process.    i K n , we can derive the following sgCME in compact form as [27]    

Generalized Gillespie algorithm
Recall that traditional GA introduces two random numbers: one for modeling which reaction will happen and the other for modeling when the next reaction will happen. The GA can provide statistically exact methods for simulating stochastic dynamics as interacting sequences of discrete events [29]. However, this algorithm exploits the hypothesis that inter-event times follow exponentially distributions. In the present study, we relax this hypothesis and propose a so-called sgGA based on Eq. (4), which is an extended version of the GA but can be applied to stochastic analysis of a general biochemical network with arbitrary waiting-time distributions. It is worth pointing out that many other extensions of the GA can also be established in cases of non-Markovian 8 reactions.
Similar to the traditional GA, we draw two random numbers 1 r and 2 r from the uniform distribution in the unit interval, and set After discarding the early part of the trajectories (the burn-in phase), the remaining realization values are used to estimate the steady-state probability distribution and related statistics.

Generalized finite state projection algorithm
The FSP method developed initially by Munsky and Khammash has been very successful in solving the CME even for relatively large systems [45]. Based on the idea of FSP, here we propose a sgFSP, which enables to compute stationary probability distributions in systems of non-Markovian biochemical networks.
First, the sgCME (i.e., Eq. 4) can be rewritten in the following form: 0  AP (6) where all possible states space are enumerated as a vector   We point out that if the number of the states is infinite or extremely large, the sgCME is very difficult to solve, but one can get a good approximation of the solution by making use of truncation methods such as finite state projection [45], sliding window [46], finite buffer method [47,48] and quantized tensor trains methods [49]. In these methods, truncation is made to guarantee both that the number of states retained is small enough and that the truncated equation is efficiently solved.
On contrary, if the number of states retained is very large, the truncated equation is still difficult to solve.
Then, we adopt finite state projection as a foundation to construct a sgFSP. Let denote a Markov chain on the state set , whose master equation is given by  10 The first example to be analyzed is a birth-death model, which is a natural and formal framework for modeling a vast variety of biological processes such as population dynamics, speciation, genome evolution [27,50]. Consider a birth-death network consisting of the following two non-Markovian reactions (Fig. 2(A))

A model of birth and death with molecular memory
where two waiting time distributions are set as i.e.,  is the common Gamma function. In numerical simulation, we limit ours consideration to the case of Gamma distributions, although our approach may be applied to other kinds of distributions suggested in the literature [51] (additional explanations can be found in the next section).
It is worth pointing out that the above formulation contains the stochastic model of a gene's constitutive expression as its special case, where protein synthesis corresponds to birth while protein degradation to death [52,53]. Here, we consider the following two cases of gene expression (referring to Figure 2(C) and (F)):  Fig. 2(B)   Numerical results are demonstrated in Fig. 2. In Fig. 2(D)(E), lines correspond to distributions obtained by sgFSP and empty circles to distributions obtained by sgGA. From these two panels, we observe that the two obtained kinds of distributions match remarkably well, implying that the two methods are equivalent. Therefore, we focus on the performance of sgFSP in the following. We also observe that protein distributions shifts toward the direction that b k increases (Fig. 2(D)(G)(H)), but toward the opposite direction ( Fig.2 (E)). Interestingly, protein populations demonstrate bimodal distribution when 1 b k  (Fig.2 (E)), implying that molecular memory can induce population diversity.

Gene model of toggle switch with molecular memory
Recall that a toggle switch network (Fig. 3(A)) can model the cross-repression between the determinants of different cellular states, which can result in a definite choice between two outcomes [54][55][56]. Conventional models of genetic toggle switch consider exponential waiting time distributions. However, the expression of a gene in general involves a multi-steps process. Indeed, transcriptional repressor monomer (A or B) binds first to dimers and then to specific DNA sequences near the promoter, repressing the production of transcriptional repressor monomer (B or A). This multi-step process can result in non-exponential waiting times and hence create molecular memories between reaction events. Here, we consider a stochastic model of genetic toggle switch with memories, which consists of the following four reactions (or referring to Fig. 3 Fig. 3 (H), we find that two switching states occur only in the case of non-exponential waiting times or memory. In a word, we conclude from Fig. 3 that molecular memory can induce bimodal distributions in the toggle switch model described in Fig. 3(A) or Fig. 3(B).

On-off model of gene expression with molecular memory
The 3-stage model of gene expression [12][13][14][15]57,58] can find its prototypes in prokaryotic and eukaryotic cells. This model, which can capture many essential features of transcriptional biology, has extensively been used and studied. The basic hypothesis of the model is that the gene promoter exists either in an activated state that produces mRNA or in a repressed state that is unproductive, and that an mRNA undergoes a competition between translation and degradation. In traditional studies, it is also assumed that every chemical reaction involved happens probabilistically at a fixed rate, i.e., waiting times between reaction events are assumed to follow exponential distributions.
However, biochemical events do not necessarily occur in single steps of individual molecules. For example, a promoter's gradual maturing through a series of inactive states (many of which have been unspecified) before activation can result in non-exponential distributed gestation periods between transcription windows. Also for example, mRNAs or proteins' senescing through a series of states (many of which have also been unspecified) before they are finally degraded can also lead to non-exponential waiting times. In order to model various possible regulations involved in gene expression, we replace exponential waiting time distributions with non-exponential waiting time distributions for transitions between promoter states, transcription, translation and degradation. Thus, the traditional 3-stage model of gene expression is modified as (also referring to Fig. 4(A)): We assume that all these waiting time distributions are Gamma distributions, that is, (1) Promoter switch process: (2) Transcription process: As is well known, the relative rate of the constant gene-state transition rates has significant impact on gene expression dynamics. Here, we compute protein distributions in two dynamic regimes: 'fast' gene-state dynamics (possibly due to fast binding and unbinding of transcription factors)) and 'slow' gene-state dynamics (possibly due to slower changes in chromatin configuration). Note that parameters on L and off L in the model represent the step numbers of multistep processes, which result in Gamma distributed waiting times. In simulation, the average waiting times is kept constant. Numerical results are shown in Fig.4 (B)-(G). From these panels, we observe that whichever "fast" and "slow" switches, protein distributions shifts to smaller number density as on L is increased (Fig.4 (B)(E)), but shifts higher number density as   off tot LL is increased ( Fig.4 (C)(D)(F)(G))). Interestingly, protein populations tend to limit distributions as

Discussion
A key challenge in stochastic modeling of cellular networks is in raising computational efficiency involving the jointing effects of stochastic fluctuations (or the noise) and molecular memory. We have developed a novel methodological approach, which is based on sgCME, to compute steady- First, our methods rest on sgCME. Although this equation can reveal long time behavior of stochastic process with memory, but lacks moderately transient dynamics. To overcome this difficulty, there are two possible roadmaps: One is that one may take the whole line back to the gCME formulation, numerical integrate the corresponding integro-differential equation, and store the whole historical data; And another is to resort to Monte Carlo simulation. Future extensions of the algorithms are under way. Second, the waiting time distribution selection is essentially a statistic inference problem in systems biology. Large amount of experimental data imply that Gamma distribution and phase-type distribution are good choices. The specifics are outside the scope of this article and will be discussed elsewhere.
We stress that our framework opens a door to analyze and simulate stochastic behaviors of the systems of cellular processes with molecular memory. The presented formulation enables us to efficiently develop algorithms under a Markovian framework with effective transition rates instead of reaction propensity functions. And our theory is general in that it can deal with biochemical networks of arbitrary topology and arbitrary types of waiting time distributions. Therefore, theoretically, successful analysis methods and simulation algorithms developed in the Markovian framework can apply to non-Markovian problems [59]. In addition, it is worth pointing out that our 18 formalism is not restricted to biochemical networks considered here and can be easily extended to any non-Markovian networks with a finite/infinite set of discrete states, such as epidemic spreading on complex networks [60][61][62][63].
Finally, the emergence of measurement techniques currently allow for access to increasingly rich data on approximately steady-state distributions of interested molecular species and waiting time distributions of key reaction events. The presented methods are computationally efficient and scalable. This will facilitate the quantitative mechanistic modelling of complex cellular processes and the exploitation of cell-to-cell variability and molecular memory.

Appendix: Biophysical foundation of waiting-time distributions
The timing of cellular processes cannot be described often by simple Poisson processes, which are traditionally characterized by exponential distributed waiting time distributions and hence can be modelled by classical chemical master equations. However, many single-peaked and even multimodal waiting-time distributions have been observed empirically, reflecting complexity of cellular processes modelled by biochemical reactions networks. Indeed, non-exponential waiting-time distributions are signatures of relevant intracellular processes without a priori knowledge. A main advantage of gCME capturing effects of non-exponential waiting times is that it can focus on a small number of key measurable events. Here, we derive the most common waiting time distributions of biochemical reactions: exponential distribution, Gamma distribution and phase-type distribution, focusing on their biophysical explanations.

Exponential distribution
Exponential distributions can characterize the time between events in a Poisson point process.
Exponential waiting times between reactions are a pillar of the classical chemical master equation. It is the exponential distribution which can be used to modeling the waiting times between single elementary reaction steps. Note that this derivation is based on Markovian assumption and the exponential distribution can also be derived from a master equation [3,4].

Gamma distribution
In contrast to the times until the next molecular event in single-step processes (elementary reactions), the waiting times between reactions are, in general, not exponentially distributed. If a reaction process goes through a series of (identical) exponential steps, the waiting times will be Gamma distributed [51]. Mathematically, this process can be regarded as a sum of n single-step processes, and therefore the over-all waiting time is the n-fold convolution     where, * denotes convolution and   t  is known as the gamma distribution with shape parameter n and rate parameter μ. If n is a positive integer, then the distribution is an Erlang distribution.

Phase-type distribution
A phase-type distribution is a probability distribution constructed by a convolution or mixture of exponential distributions. It results from a system of one or more inter-related Poisson processes occurring in sequence or phases, such as reversible chemical reactions. For example, models 20 describing enzyme dynamics can be mapped onto phase type models. Phase-type distributions can provide a tractable theoretical framework to understand complex waiting time distributions observed for a variety of different biochemical processes. Formally, the PDF of a phase-type distribution is [51]  There are a lot of probability distributions involving waiting times, such as Weibull distribution, Pareto distribution, and Mittag-Leffler distribution. We do not tend to enumerate all here.