Exact Probability Distributions of Selected Species in Stochastic Chemical Reaction Networks

Chemical reactions are discrete, stochastic events. As such, the species’ molecular numbers can be described by an associated master equation. However, handling such an equation may become difficult due to the large size of reaction networks. A commonly used approach to forecast the behaviour of reaction networks is to perform computational simulations of such systems and analyse their outcome statistically. This approach, however, might require high computational costs to provide accurate results. In this paper we opt for an analytical approach to obtain the time-dependent solution of the Chemical Master Equation for selected species in a general reaction network. When the reaction networks are composed exclusively of zeroth and first-order reactions, this analytical approach significantly alleviates the computational burden required by simulation-based methods. By building upon these analytical solutions, we analyse a general monomolecular reaction network with an arbitrary number of species to obtain the exact marginal probability distribution for selected species. Additionally, we study two particular topologies of monomolecular reaction networks, namely (i) an unbranched chain of monomolecular reactions with and without synthesis and degradation reactions and (ii) a circular chain of monomolecular reactions. We illustrate our methodology and alternative ways to use it for non-linear systems by analysing a protein autoactivation mechanism. Later, we compare the computational load required for the implementation of our results and a pure computational approach to analyse an unbranched chain of monomolecular reactions. Finally, we study calcium ions gates in the sarco/endoplasmic reticulum mediated by ryanodine receptors. Electronic supplementary material The online version of this article (doi:10.1007/s11538-014-9985-z) contains supplementary material, which is available to authorized users.


Solution of an ODE with Repeated Eigenvalues
Here we derive the closed form solution of a linear system of the form d dt e(t) = Ae(t), e(0) = e 0 , where e(t) : R + → R n , A ∈ R n×n , and C selects a number of states or linear combinations of e(t). We consider the case in which A has repeated eigenvalues.
To obtain the closed-form solution of (1), we express its transfer function as a finite sum of fractions. This transfer function is given by C (ζI − A) −1 and describes the relationship in the frequency domain between initial conditions e 0 and output y(ζ); that is to say Here y(ζ) is the Laplace transform of y(t) and ζ represents the complex frequency variable that arises from the Laplace transform. The following proposition describes how to express such a transfer function as a finite sum of fractions.
Proposition 1. (Leyva-Ramos, 1991) Consider a system of the form (1), where the matrix A has σ eigenvalues with multiplicity µ ∈ Z ∀ ∈ [1, σ]. The partial fraction expansion of C (ζI − A) −1 is given by where  ⊗ is the Kronecker product, and V is the Vandermonde matrix formed from the eigenvalues of A.
The Vandermonde Matrix of n numbers λ , where λ i is repeated µ i times is (Leyva-Ramos, 1991) In the equation above, the superindex denotes exponentiation. By using the partial fraction expansion in (3) and by taking the inverse Laplace transform of (2), we obtain the solution of (1): In Section 3.2 of the main manuscript, we are interested in deriving the solution of the following dynamical systems Here C is a matrix of appropriate dimensions that select the entries of interest of q j and ν; and ε j is the j th column of an identity matrix. When A has σ eigenvalues each with multiplicity µ , we can avail of Proposition 1 of this SOM to obtaiñ where K i is given in (4), k j i is the j th column of K i , andν represents the steady state of (5b).

Steady-State Marginal Probability Distribution for an Unbranched Monomolecular Reaction Chain
This section is a follow-up of the case-study in Section 4.2 of the main manuscript. The reaction network studied is Here, we derive steady state distribution by two different methods: (i) the analytical expressions in the (23) and Proposition 2 of the main manuscript and (ii) the multiscale methodology by Cotter et al. (2011), that assumes that the evolution of the slow species is well approximated by a Langevine process. For our analytical approach, we note that the steady-state probability density function described in Proposition 2 of the main manuscript is given by the expression lim t→∞ F (y, t) = P (y,ν(∞)) * M y, ξ 1 ,q 1 (∞) * . . . * M (y, ξ n ,q n (∞)) , whereq j (∞) andν(∞) is an abuse of notation, which we use to denote the steady state of the ODEs in (5). In turn, y gathers the molecular species number of the first and last species, i.e. y := (s 1 s n ) T . When all the kinetic parameters of (7) are positive, the eigenvalues of the matrix A in (5) are all negative. Hence from (23) of the main manuscript, we conclude that The matrix C is designed to obtain the first and last entry of the vector A −1 b. By using these expressions in the steady-state distribution above, we obtain lim t→∞ F (y, t) = P y, −CA −1 b .
We note that this is the exact solution for the steady state of the chemical master equation for the selected species. Now, the multiscale methodology in (Cotter et al., 2011) requires the identification of slow and fast varying species. Once this separation of scales is identified, the probability distributions that characterises the behaviour of the fast variables is obtained via stochastic simulations for each possible combination of values of the slow variables. Please, refer to (Cotter et al., 2011, Sec. III) for a detailed explanation. The outcome of these simulations is used to approximate the steady-state of the slow species via the steady-state of a Fokker-Planck equation. To obtain the solution of such a Fokker-Planck equation, we availed of the Finite Volume Method described in (Sjöberg et al., 2009, Sec. 2.2). Figure 1 depicts the comparison of the outcome of both approaches.
Finally, we compared the computational time required to obtain such a steady-state distribution by means of the expression Here t MS denotes the computational time required by the method in (Cotter et al., 2011) and t A represents the time required by the analytical method. Figure 2 depicts such a computational time comparison for different lengths n of the reaction network (7).